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Abstract 

The surrounding world surprises us by the beauty and variety of complex shapes 
that emerge from nanometric to macroscopic scales. Natural or manufactured ma- 
terials (sandstones, sedimentary rocks and cement), colloidal solutions (proteins 
and DNA), biological cells, tissues and organs (lungs, kidneys and placenta), they 
all present irregularly shaped "scenes" for a fundamental transport "performance" , 
that is, diffusion. Here, the geometrical complexity, entangled with the stochastic 
character of diffusive motion, results in numerous fascinating and sometimes un- 
expected effects like diffusion screening or localization. These effects control many 
diffusion-mediated processes that play an important role in heterogeneous catalysis 
(reaction rate and overall production), biochemical mechanisms (protein search and 
migration), electrochemistry (electrode-electrolyte impedance), growth phenomena 
(solidification, viscous fingering), oil recovery (structure of sedimentary rocks), or 
building industry (cement hardening) . In spite of a long and rich history of academic 
and industrial research in this field, it is striking to see how little we know about 
diffusion in complex geometries, especially the one which occurs in three dimensions. 

We present our recent results on restricted diffusion. We look into the role of 
geometrical complexity at different levels, from boundary microroughness to hier- 
archical structure and connectivity of the whole diffusion-confining domain. We 
develop a new approach which consists in combining fast random walk algorithms 
with spectral tools. The main focus is on studying diffusion in model complex ge- 
ometries (von Koch boundaries, Kitaoka acinus, etc.), as well as on developing and 
testing spectral methods. We aim at extending this knowledge and at applying 
the accomplished arsenal of theoretical and numerical tools to structures found in 
nature and industry (X-ray microtomography 3D scans of a cement paste, 2D slices 
of human skin and the placenta, etc.). 
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Preface 

The essay is based on the author's research summary that has been prepared for defend- 
ing a Habilitation degree (Habilitation a Diriger des Recherches in French). This is a 
structured overview of the author's works in different fields in which diffusion, geome- 
try and complexity are entangled. As a summary of one scientist's activity the essay 
is unavoidably biased and is not meant to provide a reader with a complete overview. 
Although the bibliography is extended, it certainly does not include all the references 
that are relevant for these research fields. In spite of this incompleteness, a reader may 
hopefully find interesting results and intriguing open problems for future research. 
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Introduction 



Irregular shapes are ubiquitous in fields as different as biology, physiology, chemistry, 
electrochemistry, engineering, and material sciences. They play a crucial role in various 
natural phenomena and condition multiple industrial processes. For instance, a human 
being needs around hundred square meters of alveolar surface in the lungs for respiration. 
Since this surface has to fit the rib cage, it is necessarily irregular. Catalysts are often 
manufactured to have a rough boundary in order to increase the overall production be- 
cause of a larger reactive surface. Oil-bearing sedimentary rocks and sandstones, cement 
and concretes, biological and textile tissues exhibit porous, often multiscale structures, 
in which transport processes take place. Either intrinsic or artificially implemented, a 
geometrical complexity essentially determines the properties and the overall functioning 
of an organ or an industrial device. Even for such a well-studied transport process as 
diffusion, an irregular geometry yields new, often unexpected phenomena. The central 
question of this work is how a geometrical irregularity influences the dynamics of a com- 
plex system and, in turn, what information on geometry can be extracted from surveying 
this dynamics. 

We focus on diffusion for many reasons. First, it provides a reliable macroscopic 
framework for describing various microscopic dynamics according to the central limit 
theorem [1]. As such, diffusion turns out to be a fundamental transport mechanism in 
biology and chemistry. Second, diffusion is governed by a simple, though remarkably 
rich, second order partial differential equation (PDE) named diffusion (or heat) equation. 
Its interpretation in terms of Brownian motion provides the whole arsenal of powerful 
probabilistic tools as well as a basis for intuitive reasoning. Third, the spectral theory 
of the underlying Laplace operator is probably the most developed branch of functional 
analysis. Although diffusion is studied for a long time, the presence of a geometrical 
irregularity may drastically alter well-established results. As a matter of fact, it is striking 
to see how little we know about diffusion in complex geometries, especially the one which 
occurs in three dimensions. This is the reason for increasing interest in diffusion in recent 
decades. 

Numerous questions and problems concerning diffusion in complex geometries can be 
formally divided in three groups: 

• forward problems which consist in finding transport properties for a known (de- 
fined) geometry of a diffusion-confining domain. In mathematical terms, this is 
equivalent to finding a solution of the related PDE in this domain. More specifically, 
we are interested in finding a way to understand, both qualitatively and quantita- 
tively, how various geometrical features of the confining domain change transport 
properties. 

• inverse problems which aim at determining an unknown geometry of a diffusion- 
confining domain from measuring the transport properties of a system. This is the 
main purpose of imaging techniques. Diffusion-weighted nuclear magnetic resonance 
(NMR) imaging is an example. In this case, the magnetization attenuation which 
is caused by restricted diffusion in a porous medium is measured in order to get its 
global geometrical characteristics like surface-to- volume ratio, pore size distribution, 
etc. Kac's question "Can one hear the shape of a drum?" is the best illustration 
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of these issues [2]. Given their practical importance (e.g., in oil recovery or medical 
diagnosis), we shall bear in mind inverse problems in the course of this work. 

• optimization problems which focus on conception and design of an optimal ge- 
ometry to maximize or minimize certain transport characteristics of a system. The 
examples are efficient heat radiators, gas exchange units (e.g., filters), or catalysts, 
to name a few. Although we are not currently studying these questions, optimization 
problems are considered as an important direction for future research. 

Throughout this work, we consider forward problems by adapting two approaches: 
Monte Carlo simulations and spectral analysis. As for Monte Carlo techniques, they 
prove to be efficient and very flexible tools for generating individual random trajectories 
of diffusing particles in complex geometries, including fractal boundaries, porous media, 
multiscale or branching structures, etc. (Section 2). As far as spectral analysis is con- 
cerned, it provides a fundamental explanation and a clearer interpretation of the observed 
effects (Section 3). The goal here is to relate various features of restricted diffusion to 
the Laplace operator eigenbasis and then to investigate the latter in depth. Since the 
numerical computation of the Laplace operator eigenbasis is a difficult task for complex 
geometries, we mainly focused on simple domains (e.g., interval, disk, sphere). For these, 
many characteristics of restricted diffusion can be found either analytically, or very ac- 
curately via spectral-oriented numerical techniques. As for complex geometries, their 
systematic study is one of the principal research lines in future. 

The manuscript is organized as follows. Section 1 briefly summarizes physical concepts 
and mathematical tools for describing diffusion in complex media. We defined two prin- 
cipal research guidelines in order to present our main results: Monte Carlo simulations 
(Section 2) and spectral analysis (Section 3). All the author's publications are available 
online 1 . 



^ttp : //pmc . polytechnique . f r/pagesperso/dg/publi/publi_e .htm 
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1 Diffusion, geometry and complexity 



1.1 Diffusion 

Diffusion is a fundamental transport mechanism, with countless examples in nature and 
applications in sciences, from physics to biology, chemistry, medicine, engineering, and 
economics. 

"Upscaling" : toward macroscopic (or PDE) description 

The microscopic dynamics of gases and liquids can be described by kinetic theory, when 
a collision integral accounts for short-time and short-range interactions between atoms or 
molecules [3, 4]. From this classical point of view, the trajectory of a particle is a sequence 
of "moves" between successive "collisions" (or effective interactions) with other particles 
or a reservoir. In many cases, the mean free path is much shorter than macroscopic length 
scales at which the system is being observed. In addition, a large number of interactions 
considerably reduces or fully eliminates memory effects at macroscopic time scales, so 
that the specific motional features of the particle trajectory are averaged out. Performing 
such an "upscaling" , one aims to introduce a macroscopic density of particles and study 
its properties instead of looking at the microscopic dynamics of these particles. While 
the microscopic dynamics may be defiantly complex, its coarser macroscopic description 
is often appropriate and very accurate. 

Diffusion equation 

The diffusion (or heat) equation provides the simplest, though very general, coarser de- 
scription of the microscopic dynamics: 

^c(T,t)=DAc(r,t), (1) 

where A = d 2 jdx\ + ... + d 2 /dx 2 l is the Laplace operator in d dimensions, and D the free 
diffusion coefficient. This equation states that the evolution of the density of particles 
c(r, t) in time (the left-hand side term) is only caused by local displacements of particles 
in space (the right-hand side term). The differential form of this equation reflects the local 
character, in space and time, of the microscopic dynamics. If large moves were allowed 
during a short time, the Laplace operator A would be replaced by an integral operator 
accounting for particles coming from distant regions (spatially nonlocal macroscopic dy- 
namics). If the microscopic moves were highly correlated, the time derivative d/dt would 
be replaced by an integral operator accounting for the memory of the whole preceding 
evolution (temporarily nonlocal macroscopic dynamics). Although both situations can be 
encountered, we shall not consider such anomalous diffusions (see [5-10] and references 
therein). 

Diffusive propagator 

The diffusion equation describes the time evolution of the spatial density c(r, t) from a 
given initial state c(r = r ,t = 0) = p(r ). Among various initial conditions, a point-like 
source plays a special role. In fact, it is natural to ask how the density of particles started 
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from a given point r evolves? The family G t (r ,r) of these densities, parameterized by 
the starting point r , bears different names: diffusive propagator, heat kernel, or Green 
function of diffusion equation. These densities satisfy the diffusion equation with respect 
to r, 



with the initial condition for all particles to be concentrated in one point r which is 
mathematically represented by the Dirac distribution 



The symmetry property Gt(r, r') = Gt{r',r) holds. 

The propagator is an elementary "block" describing the dynamics of particles at macro- 
scopic level. As such, the propagator contains all the available information about diffusion. 
In particular, the solution c(r, i) of diffusion equation (1) with a given initial density p(r ) 
is simply expressed through the propagator: 



Extensions 

Apart from the aforementioned global modifications (spatially or temporarily nonlocal 
dynamics), the diffusion equation (1) can be extended in many ways: 

1. In anisotropic media, diffusion in some spatial directions may be faster or slower 
than in the others. This can be easily incorporated by adding weighting coefficients 
in front of the partial derivatives d 2 jdx\. In general, the Laplace operator can be 
replaced by a second-order elliptic differential operator, 



where the coefficients aij(r) satisfy some requirements to ensure the ellipticity [11]. 
Clearly, this extension includes a spatially inhomogeneous diffusion coefficient. It 
is worth noting that many results for the Laplace operator relying on the spectral 
analysis are directly applicable to this extension. In what follows, we shall focus on 
the Laplace operator. 

2. In some cases, the diffusion coefficient depends on time t or/and density c(r, t). The 
latter dependence may describe saturation effects in catalytic reactions. However, 
this dependence would make the diffusion equation nonlinear, significantly reducing 
the arsenal of analytical tools. For this reason, we shall not consider such extensions. 

3. A bulk source or sink term can be added to the diffusion equation. Its practical 
implementation is straightforward, especially when using the spectral description. 

4. An external scalar field B(r,t) can be included 




(2) 




(3) 





d_ 

dt 



c(r, t) = DAc(r, t) - kB(t, t)c(r, t), 



(4) 
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where k is a parameter. The last term accounts for the local effects: loss or pro- 
duction of the particles, modification of their state, etc. The proportionality to the 
density c(r, t) reflects that these local effects act on each particle individually and 
independently from the other particles. The second and higher-order powers of the 
density may represent various chemical reactions, but we shall not consider these 
nonlinear equations. 

Depending on the application field, B(r,t) incorporates various mechanisms. If 
the bulk contains absorbing sinks or traps for diffusing particles, B(r,t) represents 
the distribution of their absorption or trapping rates. If the particles possess some 
activity function which can relax on bulk impurities (e.g., a magnetization of a 
polarized spin carried by a diffusing nucleus), B(r,t) describes bulk relaxation rate. 
If the particles can be chemically transformed by catalytic germs (or biological 
cells) distributed in the bulk, B(r, t) characterizes the distribution of reaction rates. 
In NMR, B(r,t) with imaginary parameter k represents the effect of the applied 
magnetic field onto the transverse magnetization (Sect. 1.5). In what follows, our 
main focus will be onto Eq. (4). 

"Downscaling" : toward probabilistic description 

In a coarser description via diffusion equation, all the information about individual tra- 
jectories of the particles is lost. Nevertheless, thinking of diffusion in terms of individual 
particles can still be instructive to better understand or interpret various characteristics 
of diffusion. As we already mentioned earlier, the real microscopic dynamics is in general 
very complicated and intractable analytically. To overcome this problem, one employs the 
concept of "downscaling" which consists in finding a simpler microscopic dynamics that 
would lead to the chosen macroscopic description. A Brownian motion, mathematically 
defined as a stochastic process with independent Gaussian increments [1], is a very natural 
candidate for substituting the real microscopic dynamics. In fact, the diffusive propaga- 
tor Gf(ro,r), satisfying Eqs. (2, 3), turns out to be the probability density for Brownian 
motion to move from r to r during time t. This is a direct link between probabilistic and 
PDE descriptions. Strictly speaking, this "substitution" is neither unique, nor physically 
justified. It is rather a useful model that offers efficient probabilistic tools for theoretical 
and numerical analysis. In particular, random walks on a lattice (a discretized version of 
Brownian motion) are broadly employed to model many physical processes [12, 13]. 

Feynman-Kac formula 

In practice, it is difficult to monitor the individual Brownian trajectories of the diffusing 
particles. It is much easier to measure or investigate "an observable", or a functional, of 
these trajectories. We consider the random variable 



where Brownian motion is started with a given initial density p(r ). Intuitively, a given 
function B(r,t) can be thought of as a distribution of "markers" to "color" the trajectory 
X t exploring different points or specific regions in space. When the diffusing species 




(5) 



o 
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passes through these regions, the random variable tp accumulates the corresponding marks. 
In other words, different parts of the trajectory are weighted according to the function 
B(r, t), encoding thus the whole stochastic process. For example, if B(r, t) represents the 
distribution of absorption or relaxation rates in the bulk, tp is the cumulant absorption 
factor penalizing the trajectories that pass through the sinks or traps. In NMR, the 
encoding mechanism is experimentally realized by applying an inhomogeneous magnetic 
field B(r,t), and tp is the total phase accumulated by an individual spin-bearing particle 
during its diffusion in this field [14]. 

Since tp is a random variable, it is convenient to consider the Fourier or Laplace trans- 
form of its probability distribution which can be respectively interpreted as characteristic 
function Eje* 13 ^} and survival probability K{e~ Kip }. The Feynman-Kac formula relates 
these probabilistic quantities to the solution of Eq. (4) 



(if k is replaced by —iq, one obtains the characteristic function) [15-19]. In what fol- 
lows, we shall frequently employ this relation to "switch" between probabilistic and PDE 
representations. In particular, the probabilistic description suggests using Monte Carlo 
simulations for solving PDE problems [20-22]. We shall consider various applications of 
this numerical technique in Section 2. 

1.2 Boundary condition 

When the motion of diffusing species is restricted inside a confining medium, physico- 
chemical or biological interactions between the particles and the interface of the medium 
should be taken into account. For instance, paramagnetic impurities dispersed on the 
interface cause surface relaxation in NMR experiments; cellular membranes allow for a 
semi-permeable transport through the boundary; chemical reaction may change the dif- 
fusive or magnetic properties of the particle, etc. A reliable description of these processes 
at microscopic level is a challenging problem, demanding for example accurate molecular 
dynamics simulations near the interface, or quantum mechanics calculations. At the time 
scale of the macroscopic transport process, however, the contact with the interface is very 
rapid so that the precise description of the interaction is often irrelevant. In analogy with 
diffusion coefficient D, effectively representing the bulk dynamics, the interactions on the 
boundary can macroscopically be described by a surface transport coefficient W [23, 24]. 
In microbiology, this is the permeability characterizing the rate of transfer across a semi- 
permeable membrane. In heterogeneous catalysis, W is the reactivity of a catalyst, that 
is the rate at which diffusing species are chemically transformed into other species after 
hitting the catalytic surface. In NMR, W is the surface relaxivity determining the rate 
at which the nuclei lose their magnetization in the vicinity of the boundary. In the two 
latter cases, the transformed or relaxed species still remain inside the confining domain 
but they do not participate at the transport process any more (e.g., they do not contribute 
to formation of the macroscopic signal, see Sect. 1.5). 

At macroscopic level, interactions with the interface are effectively incorporated via a 
boundary condition, which is a mass conservation law: the flux of particles from the bulk 
towards the boundary, —Ddc(r,t)/dn, is equal to the flux through the boundary (or the 




(6) 
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flux of transformed or relaxed species) Wc(r,t), yielding Robin (also known as Fourier, 
mixed, relaxing, radiative, or third) boundary condition 



D—c(T,t) + Wc(r,t) = 0, 



(7) 



where d/dn is the normal derivative pointing outward the bulk (if the domain is a sphere, 
the normal derivative is equal to the radial derivative). When the surface transport 
coefficient W is zero (no flux through the boundary), one retrieves Neumann boundary 
condition: dc(r,t)/dn = 0. The opposite limit of infinite W (no resistance to the transfer 
through the boundary) yields Dirichlet boundary condition: c(r, t) = 0. The intermediate 
Robin boundary condition is a linear combination of these two extreme cases, which are 
"weighted" by bulk and surface transport coefficients D and W , respectively. In other 
words, the ratio D/W represents "proportions" of pure reflections (the first term) and 
pure absorptions (the second term), mixing Neumann and Dirichlet boundary conditions. 

The role of Robin boundary condition for Laplacian transport was thoroughly in- 
vestigated (see [23, 24] and references therein). In particular, the ratio D/W , which is 
homogeneous to a length, bears different names: "unscreened perimeter length" or "ex- 
ploration perimeter or length" [25], "Damkohler first ratio" in chemistry [26], "surface 
relaxation length" in NMR [27]. In the latter case, D/W is the distance a particle should 
travel near the boundary before surface relaxation reduces its expected magnetization. In 
Sect. 2.3, we shall give another interpretation to this exploration perimeter. 

1.3 Reflected Brownian motion 

In the above macroscopic description, restricting walls of a confining medium were intro- 
duced through boundary conditions for diffusion equation. Alternatively, this effect can 
be incorporated for Brownian motion in the probabilistic description. We first consider a 
much simpler situation of fully absorbing interface (Dirichlet boundary condition). For a 
Brownian motion W t started from a point r inside a confining domain Q, we define the 
first hitting time T when W t reaches the boundary dVL of Vt: 



Since the interface is fully absorbing, the process is stopped (or killed) at this (random) 
moment. Figuratively speaking, we simply "close our eyes" on what happens with Brow- 
nian motion after T. It is relatively easy to incorporate Dirichlet boundary condition by 
restricting the analysis to times t inferior to T. 

The situation with a reflecting interface (Neumann boundary condition) is completely 
different. In this case, one needs to keep Brownian motion inside the confining domain 
without stopping it. In other words, we have to modify the local dynamics of this process 
to include reflections on the boundary. In sharp contrast with ordinary Brownian motion, 
the construction of reflected Brownian motion strongly depends on the geometry of a 
confining medium, being especially sophisticated for irregular boundaries. The simplest 
example of this process is the absolute value of a one-dimensional Brownian motion Wf. 
X t = \Wt\. Each time a diffusing particle crosses the boundary of the positive semi- 
axis (the end point 0), it is reflected toward the bulk (positive semi-axis). In a general 
situation of a bounded domain Q C M. d with twice continuously differentiate boundary 



T = w£{t > : W t G dtt}. 



(8) 
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dfl, reflected Brownian motion X t can be denned as a solution of the stochastic differential 
equation, called Skorokhod equation [17]: 

dX t = dW t + n(X t )I dn (X t )de t , 

where W t is the (ordinary) Brownian motion, n(r) the normal inward vector at boundary 
point r, lan(r) the indicator function of the boundary (Ign(r) = 1 if r e dQ, and 
otherwise), and It the local boundary time process, satisfying certain conditions [17]. The 
most unusual feature of this definition is that the single equation defines two processes, 
X t and £ t , strongly dependent on each other. The intuitive meaning of the Skorokhod 
equation is simple. In the bulk, an infinitesimal variation dX t of the reflected Brownian 
motion inside the confining domain is only governed by the variation dWt of the (ordinary) 
Brownian motion (the second term vanishes due to the indicator function Ion)- When the 
particle hits the boundary, the second term does not allow to leave the domain leading to 
a variation directed along the inward unit normal n(r) towards the interior of the domain. 
At the same time, each encounter with the boundary increases the local time £ t . Being 
defined in this way, reflected Brownian motion can be related to diffusion equation with 
Neumann boundary condition. It is worth noting that reflected Brownian motion can 
alternatively be introduced as a continuous limit of reflected random walks on a regular 
lattice, which are easier for intuitive interpretation [28, 29]. We also mention that reflected 
Brownian motion can be rigorously defined for domains with irregular boundaries using 
the Dirichlet form method [30, 31]. 

In the intermediate case of partially absorbing/reflecting interface (Robin boundary 
condition), one can introduce "partially reflected Brownian motion" [32]. This is reflected 
Brownian motion which is conditioned to stop (i.e., to be absorbed) on the boundary at 
random moment T h when the local time £ t exceeds an independent exponentially dis- 
tributed random variable £: 

T h = inf{t>0 : £ t > £}, where P{£ > A} = exp[-A/i]. 

The absorption can happen whenever partially reflected Brownian motion hits the bound- 
ary, these hits being counted by the local time. A positive parameter h, a kind of ab- 
sorption rate, is actually the ratio between the surface and bulk transport coefficients, 
h = LW/D, L being a characteristic size of the domain to get rid off the dimensional 
units. In turn, the exponential character of the random "threshold" £ is related to the 
fact that the absorption event is independent from hit to hit (as the exponential decay of 
radioactive nuclei). 

When h goes to infinity (Dirichlet boundary condition, W = oo), one gets £ = 
with probability 1, so that the stopping time describes the first moment when the 
local time £ t exceeds 0. This is exactly the moment when partially reflected Brownian 
motion hits the boundary for the first time. One thus retrieves the above definition (8) 
for a purely absorbing boundary. In the opposite limit h — > 0, £ = oo with probability 
1, yielding Neumann boundary condition (W = 0) as expected. Varying h allows one 
to explore various situations between pure absorptions and pure reflections that makes 
partially reflected Brownian motion to be a rich and flexible microscopic model for diffusive 
transport. 
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1.4 Spectral analysis 

To the diffusion (or heat) equation is associated a set of solutions for a given confining 
domain. A particular solution is fixed by setting the initial condition. On the other hand, 
one can remain in a general frame and study the diffusion equation or, equivalently, the 
Laplace operator. In this so-called spectral approach, one is searching for the Laplacian 
eigenfunctions and eigenvalues. Instead of looking for one particular solution, one aims to 
find all the "relevant" solutions at once. It is not therefore surprising that the computation 
of the eigenbasis, fully describing the operator, is a difficult task. But once it is worked 
out, any solution of diffusion equation and any property related to the Laplace operator 
can be deduced. 

The eigenfunctions u m (r) and eigenvalues \ m of the Laplace operator A in a bounded 
domain Q are defined as 

A« m (r) + A m « m (r) = (rGfl), (9a) 

D^-u m (r) + Wu m (r) = (r e dQ), (9b) 
on 

with an integer index m — 0, 1, 2, .... The eigenvalues \ m are nonnegative, the eigenfunc- 
tions u m (r) are orthogonal: 

J dr u m (r) u* m ,(r) = 5 m>m >, (10) 
n 

where the asterisk denotes complex conjugate, and 5 m , m i is the Kronecker symbol (<5 m , m ' = 
1 for m = m', and otherwise). Note that we explicitly fixed the normalization in Eq. (10). 
The diffusive propagator has the following spectral decomposition [33-35] 

G t (r,r')=^2<(r)u m (r')e' D ^. (11) 

m 

One can easily check that this sum satisfies the diffusion equation (2), while the initial 
condition (3) at t = is guaranteed by the completeness relation 

(J(r-r') = J2 u m(r) 

m 

At last, the boundary condition for the diffusive propagator follows immediately from 
Eq. (9b). 

It is clear that eigenfunctions and eigenvalues explicitly determine the diffusive prop- 
agator. The opposite is also true: the propagator can be formally used to define the 
eigenfunctions and eigenvalues. In fact, the same amount of information about diffusive 
motion is "stored" differently in the eigenbasis and in the propagator, the latter mixing 
this information in a specific way. 

1.5 Applications in NMR 

Observation of translational dynamics requires a kind of "marking" or "labeling" of the 
traveling particles for tracking their displacements in space. A magnetic field is a superb 
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experimental tool for encoding the motion of spin-bearing particles. For nuclei with 
spin 1/2 (e.g., protons of water), the magnetic field induces energy splitting into two 
levels. At thermal equilibrium, the population of the nuclei at the lower energy (i.e., 
the spins parallel with the magnetic field) is bigger than the population of the nuclei at 
the higher energy (i.e., the spins antiparallel with the magnetic field). The difference in 
these populations creates local magnetization which is parallel with the magnetic field 
(conventially, the z axis). This magnetization can be flipped into the transverse xy plane 
by a 90° radio-frequency (rf) magnetic field pulse (Fig. 1). The spins at position r start 
to precess with the Larmor frequency 7-B(r, t) which is proportional to the magnetic field 
B(r,t), 7 being the gyromagnetic ratio (a fundamental constant of the nucleus). The use 
of a spatially inhomogeneous magnetic field allows one to distinguish points or regions in 
space, as the nuclei precess faster or slower in different regions. This mechanism is widely 
applied in experiments to monitor the translational dynamics and to access the geometry 
of a confining medium [36]. 

The x and y projections of the transverse magnetization are conventionally treated as 
real and imaginary parts of a complex- valued "density" c(r, t). The magnetic-field encod- 
ing is then incorporated through Eq. (4) with an imaginary parameter k = ry, the last 
term in this equation being responsible for describing precession of the transverse mag- 
netization. This equation, known as Bloch-Torrey equation [37], describes the evolution 
of c(r, t) according to two independent mechanisms 

• diffusive migration of the spin-bearing particles (from r to neighboring points), rep- 
resented by the Laplace operator and characterized by the free diffusion coefficient 



• magnetic field encoding, when the spins at r acquire the phase shift 7i?(r,t)r, 
resulting from their precession in the transverse plane during a short time r. 

The macroscopic signal at time t is formed by the whole ensemble of the spins diffusing 
inside the confining domain Q: 



where p(r) is a sampling or pickup function of the measuring coil or antenna (usually one 
tries to get p(r) as uniform as possible). In analogy to Feynman-Kac formula (6), the signal 
can alternatively be written as the characteristic function E{e l7</3 } of the random dephasing 
ip of an individual nucleus diffusing in an applied magnetic field B{r,i) (the sampling 
function p(r) being implicitly incorporated here in the expectation). Both formulations 
are widely used in NMR literature (see [14] and references therein). 

In most practical situations, the encoding term B(r, t) is a superposition B + f(t)(g-r) 
of a static magnetic field B and a linear magnetic field gradient g, whose dependence 
on time is represented through a dimensionless temporal profile f(t). The form of the 
temporal profile can be easily varied in modern MR scanners. The simplest choice f(t) = 1 
corresponds to a free induction decay (FID) in a constant gradient. Two identical gradient 
pulses of opposite directions can be applied to form a gradient echo. If the nuclei were 
immobile, their dephasing by the first gradient pulse would be fully compensated by 
the second gradient pulse. When the nuclei diffuse, they experience various magnetic 
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Figure 1: Schematic illustration of echo formation. In a static field B along the z axis, the 
spins precess around this axis with the Larmor frequency ujq = 7-Bo> and its magnetization 
m is directed along z. (a) At time t = 0, one applies a periodic magnetic field Bi rotating 
in the transverse plane xy with frequency u (so-called n/2 or 90° radio-frequency (rf) 
pulse), (b) Still precessing, the magnetization is turning, linearly in time, toward the 
transverse plane, (c) At t — tv/2 (here r n = tt/^Bi)), the magnetization lies in the 
transverse plane, and the periodic magnetic field B\ ceases. During the following period 
AT = T/2 — 7V, the magnetization slowly relaxes back to the longitudinal direction 
(axis z). (d) At t — T/2 — r n /2, another magnetic field pulse is applied during the 
time (so-called rr or 180° rf pulse), (e) This pulse inverts the longitudinal direction 
of the magnetization, (f) After the 180° rf pulse being turned off at t = T/2 — r n /2, 
the slow relaxation during the subsequent time period AT returns the magnetization 
into the transverse plane. Moreover, the magnetizations of various spins are refocused 
thus forming a macroscopic signal (an echo) at t — T — r n /2. The rf pulses are very 
short allowing their durations (t w /2 and t w ) to be neglected. When the magnetic field 
is spatially inhomogeneous, the refocusing is not complete because the spins, precessing 
with varying Larmor frequencies, acquire different phase shifts. The resulting macroscopic 
signal depends on translational dynamics of spins, allowing one to survey this dynamics 
in an experiment. 
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fields and dephase differently. As a consequence, the rephasing is not complete, and the 
gradient-echo amplitude attenuation characterizes restricted diffusion (see below). 

Measuring the macroscopic signal E as a function of the experimentally controlled 
parameters (gradient, diffusion time, temporal profile, etc.), one aims to retrieve as much 
information on the dynamics and the geometry of confinement as possible. Solving this 
inverse problem necessarily requires knowing how a given geometry influences the macro- 
scopic signal. This forward problem will be considered throughout the manuscript. 

1.6 Complexity and irregular geometry 

The diffusion (or heat) equation and the underlying spectral analysis of the Laplace op- 
erator in a bounded domain is a classical problem in mathematics [34, 35]. What is 
new in our consideration of this old problem is an irregular geometry of a medium or 
its boundary. As a matter of fact, porous materials, concentrated colloidal suspensions, 
and physiological organs (such as lungs or kidneys) are examples of systems developing 
large specific surfaces. They all present a rich variety of shapes and exhibit complex mor- 
phologies on a wide range of length scales. In these systems, the interfacial confinement 
strongly influences the diffusive dynamics of Brownian particles. 

The human lung is a striking example of a complex transport system, in which the 
irregular geometry plays a central role. Its branching structure guarantees a rapid access 
of a large quantity of fresh air towards thirty thousands of pulmonary acini [38]. Each 
acinus is a porous, dichotomously branched gas exchange unit, in which oxygen molecules 
diffuse toward the alveolar membranes for further transfer to the blood. The complexity 
of this geometry was shown to lead to diffusion screening, one of the physical mechanisms 
for regulating the physiological efficiency of the lungs [25]. Many other physiological 
organs (such as kidneys, intestine, placenta) have also hierarchical multiscale structure. 
Moreover, biological systems often present geometrical complexity even at microscopic 
scales (like the surface of cells [39]). 

In material sciences, multiscale porous structures are as well ubiquitous. Typical 
examples are: sedimentary rocks, clays, sandstones, plasters, cement, etc. For instance, 
the geometrical structure of a cement paste is responsible for mechanical strength and the 
overall resistance of buildings, while its alteration in time by diffusion-mediated reactions 
is a major reason for aging, cracks, failures or even collapses [40]. A large part of worldwide 
reserve of crude oil is "imprisoned" in sedimentary rocks, and their porous morphology 
is a key factor for oil recovery [41]. A microroughness of metallic electrodes (scratches, 
abrasions, etc.) may substantially alter their transport properties and functioning [42]. 
Finally, the tortuous surface of proteins and DNA molecules, when looked at microscopic 
scales, is relevant for diffusion-mediated biochemical reactions. At a larger scale, industrial 
catalysts are made of very irregular shape to increase the catalytic surface aiming to 
improve the overall production rate [43]. 

In summary, irregular shapes are encountered in nature and sciences more often than 
one usually expects. These irregularities span a wide range of length scales (from nanome- 
ters to hundreds of meters) and exist in a variety of systems, from biology to mineral sci- 
ences. The diffusive motion of particles inside such a medium is substantially influenced 
by its complex geometry. This is a major research topic of this manuscript. 



15 



2 Probabilistic insight: Monte Carlo simulations 



In physics, speaking about diffusion in complex geometries invokes numerical analysis by 
default. In fact, only for few simple shapes (such as a disk or a sphere), diffusion equation 
possesses explicit solutions. All other domains require solving diffusion equation on a 
computer. A variety of numerical methods can be roughly divided into two groups. 

In the first group (finite differences, finite elements, boundary elements, etc.), a do- 
main and/or its boundary is discretized with a regular or adaptive mesh. The original 
continuous problem is then replaced by a set of linear equations to be solved numerically. 
The solution c(r, t) is obtained at all mesh nodes at successive time moments. Since the 
accuracy and efficiency of these deterministic numerical schemes significantly rely on the 
discretization, mesh construction for complex geometries turns out to be the key issue 
and often a limiting factor, especially in three dimensions. 

In the second group (Monte Carlo simulations), a probabilistic interpretation of PDE is 
employed. Solving the original continuous problem is replaced by modeling random trajec- 
tories of the underlying diffusion process. The solution is then obtained via Feynman-Kac 
formulas (Eq. (6) or similar) by averaging over random trajectories. Since there is no dis- 
cretization neither of the domain, nor of boundary conditions (the most subtle stage for 
the first group), Monte Carlo techniques are in general flexible and easy to implement. 
Two drawbacks should however be mentioned: a slow convergence to the solution (error 
typically decreases as l/y/N, N being the number of the simulated trajectories) and the 
local character of simulations (the solution is obtained at one spatial point). For this latter 
reason, Monte Carlo simulations are not well appropriate for getting the whole solution 
c(r,t). In contrast, this drawback is removed when one is interested in some average of 
the solution like, for instance, in Eq. (12). As illustrated in this Section, Monte Carlo 
techniques become particularly well suited, especially for studying diffusion in complex 
geometries. 

Monte Carlo simulations 

In the simplest case, one fixes a small time step r and generates a stochastic trajectory 
X t from a chosen starting point x by adding successively random increments £ n : 

X = x , AT( n+1)r = X nT + £ n (n> 0). 

For Brownian motion, all the increments are independent identically distributed random 
variables with normal (Gaussian) law, with mean zero and standard deviation a = \ / 2Dt 
(an extension to the multidimensional case is straightforward). A discretized version of 
this process is obtained by taking £ n = ±a with randomly chosen sign. These are random 
walks on a lattice with mesh a. Restricting walls can be taken into account either by 
stopping the process after their hit (Dirichlet boundary condition), or by reflecting the 
particle inside the domain (Neumann boundary condition), or by a combination of both 
(Robin boundary condition). 

Given the simplicity of this procedure, its implementation is easy for various domains, 
but its practical applications are in general limited to simple geometries. In fact, an ac- 
curate modeling of the trajectory requires to keep a typical step a (or the lattice mesh) 
smaller than the smallest geometrical feature of a confining domain. Since complex ge- 
ometries (e.g., fractals) are often multiscale, the distances that should be explored by 
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Figure 2: Illustration of a fast random walk algorithm. From the current position, one 
determines the distance to the boundary (here, for a flat horizontal segment, it is the 
height of the current point above the boundary) and draws the circle (or a sphere in 
3D). One then chooses a uniformly distributed random point on this circle and makes 
the jump. From this new position, one repeats the above steps, generating a sequence 
of points r , ri, r 2 ,..., until the particle approaches the boundary closer than a chosen 
threshold e. At this point, the simulation stops and the particle is considered as having 
hit the boundary. As illustrated on the right, each jump replaces a Brownian path inside 
the circle that speeds up simulations considerably. 

diffusing particles inside the confining domain are by orders of magnitude longer than its 
smallest geometrical features, requiring very long and time-consuming simulations for a 
single trajectory. 

To overcome this difficulty, Muller proposed to replace Brownian motion by an equiv- 
alent "spherical process" [44]. This method, often called "fast random walks", was em- 
ployed by many authors (see, e.g., [45-47]). The idea is to explore a confining domain as 
fast as possible, without violating probabilistic properties of Brownian motion. Starting 
from a given point, a diffusing particle executes a series of random jumps in the bulk. 
The jump length at each step is taken to be the distance £ between the current position 
and the boundary (Fig. 2). Since the sphere (or the disk in 2d) of radius £ centered at the 
current position resides in the interior of the confining domain, Brownian motion inside 
this sphere is not altered by the presence of the restricting walls. Figuratively speaking, 
a diffusing particle can "see" only a close neighborhood of its current position so that it 
simply does not "know" about the walls far away. The continuity of Brownian motion 
implies that the diffusing particle must intersect somewhere the frontier of this sphere 
before hitting the boundary. The rotational symmetry yields the uniform distribution of 
intersection points. It means that a lengthy simulation of Brownian trajectory inside the 
sphere can be replaced by a single random jump from its center to a uniformly distributed 
point on its frontier (Fig. 2). This "trick" drastically speeds up Monte Carlo simulations 
since at each step the largest possible exploration is performed. 

In practice, the computation is reduced to finding the distance from any interior point 
to the boundary. Depending on the studied geometry, this problem can be solved in 
different ways. In the next subsection, a geometry-adapted fast random walks (GAFRW) 
algorithm is presented. 
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Figure 3: Three generations of the quadratic von Koch curve of fractal dimension D = 
In 5/ In 3 1.465. At each iteration, one replaces all linear segments by the rescaled 
generator (first generation). 

2.1 Geometry-adapted fast random walks (GAFRW) 

Fractals are often used as a paradigm of complex domains [48-51]. On one hand, frac- 
tals are very irregular shapes which exhibit geometrical details at various length scales, 
"exploding" the classical notions of length, surface area or volume. On the other hand, 
self-similar or self-afnne hierarchical structures help to perform accurate theoretical and 
numerical analysis on fractals. In particular, the self-similarity of von Koch curves and 
surfaces allowed us for a rapid computation of the distance from any interior point to 
these boundaries [52]. 

To illustrate these ideas, we consider the quadratic von Koch curve (Fig. 3). When a 
random walker is far from the boundary, it does not "distinguish" its geometrical details. 
One can thus estimate the distance by considering the coarsest generation (Fig. 4). Get- 
ting closer and closer to the boundary, the random walker starts to recognize smaller and 
smaller geometrical details. But at the same time, when small details appear in view, the 
rest of the boundary becomes "invisible" . Consequently, one can explicitly determine the 
distance by examining only the local geometrical environment (see [52] for details). 

When the particle approaches the boundary closer than a chosen threshold, we say 
that it hits the boundary. Depending on the boundary condition and the problem at hand 
(see below), the particle is either absorbed (simulation is stopped), or reflected (simulation 
is resumed from a nearby bulk point). 

An advantage and eventual drawback of this GAFRW algorithm is the need for a 
specific implementation for each studied geometry. The algorithm was developed for the 
above quadratic von Koch curve and the cubic von Koch surface of fractal dimension 
D = In 13/ In 3, as well as for triangular von Koch curves of variable fractal dimen- 
sion (Sect. 2.2 A). In the next subsections, we shall discuss various applications of this 
algorithm. 

2.2 Multifractal properties of the harmonic measure 

The initial purpose for developing GAFRW algorithm was for studying the accessibility of 
a boundary to diffusing particles, which is of primary importance in growth and transport 
phenomena. Actually, the transfer or reactive capacity of an interface is crucially limited 
by its accessibility for Brownian motion, the effect known as diffusion screening (similar to 
electric screening in electrostatics). A particle diffusing towards an irregular interface has 
very little chance to reach the bottom of a deep "fjord" before hitting more prominent 
(and easier accessible) points. The resulting distribution of diffusive fluxes or arrival 
probabilities on the boundary is therefore very uneven, a small fraction of the boundary 
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Figure 4: Basic arrow-like cell which is divided into the rotated square and five small 
triangles (a). Once a Brownian particle arrived into the rotated square, the distance 
between its current position (full circle) and the boundary of the arrow-like cell (the 
generator) can be computed explicitly. Random jumps inside the rotated square can 
therefore be executed, until the Brownian particle either exits from the arrow-like cell, or 
enters into a small triangle. In the latter case, the Brownian particle starts to "see" the 
geometrical details of the next generation. The rescaled arrow-like cell (b) is then used 
to compute the distance to the boundary. Note that a distant source is placed on the top 
of the figure, while two vertical dotted lines delimit the interior of the confining domain. 
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Figure 5: A given boundary (here, the first generation of the quadratic von Koch curve) 
is covered by disks of diameter 5. 



receiving the overwhelming majority of the diffusing particles. The diffusion screening 
limits the overall production of species in the diffusion-limited regime of heterogeneous 
catalysis but allows one to design long- working catalysts [43]. 

Mathematically, this accessibility is characterized by a harmonic measure, which is 
defined for each (Borel) subset of the boundary as the probability for Brownian motion 
to reach the boundary of the confining domain for the first time at this subset [53]. This 
measure governs random growth processes (e.g., diffusion-limited aggregation, dendric 
growth, solidification, viscous fingering, morphogenesis), primary current distribution in 
electrochemistry, distribution of particle flows on a membrane, electric charge distribution 
on a metallic surface, etc. 

2.2.1 Multifractal analysis 

For a smooth boundary, the harmonic measure is fully characterized by its density. When 
the boundary is irregular (e.g., fractal), there is no density, though the harmonic measure 
is still well defined. In this case, one uses the multifractal dimensions D q of the harmonic 
measure to characterize its variations over different length scales. 

Let a finite generation g of a fractal boundary be covered by disjoint compact sets of 
diameter 5 (e.g., disks, spheres, cubes), as illustrated in Fig. 5. The harmonic measure at 
the scale 6 is represented by a finite set of hitting probabilities Pk,s,g- The variation of the 
harmonic measure with the scale S can be characterized by the behavior of the moments 



When 5 goes to 0, ((q,5,g) scales as 5 <yq ~ lS}Dq . For a smooth d- dimensional boundary, all 
the D q are trivially equal to d, independently of the real parameter q. Geometrical irreg- 
ularities may lead to different scaling of the harmonic measure for various q. Except for 
few classes of fractals [54-56], the multifractal dimensions of the harmonic measure need 
to be determined numerically. The GAFRW algorithm was used to simulate Brownian 
trajectories from a distant source toward von Koch boundaries and to calculate approxi- 
mately the distribution of the hitting probabilities Pk,s,g at a given scale 5. Defining the 
scale-dependent multifractal exponents as 



one can obtain the multifractal exponents D q by interpolation as g — > oo and 5^0. 

The efficiency of the GAFRW algorithm allowed us to simulate up to 10 10 trajectories, 
while the generation order could be increased up to g = 10 in 2D, and g = 7 in 3D. It 




(13) 



A.- 




jgXCgjj; a) 
(q- l)ln<5' 



20 



Dn « 1.129 



Da m 1.262 




D « 1.699 





Figure 6: The third generation of three triangular von Koch curves. 



is worth noting that these "limiting values" can be further extended by using parallel 
computation. 

2.2.2 Interpolation of the scale-dependent exponents 

First, we investigated the scale-dependent multifractal dimensions D q g >g and its interpo- 
lation as the generation order g goes to infinity and the scale 5 goes to 0. This is a priori a 
difficult problem because numerical simulations can only be realized for small generation 
orders and limited range of scales. We discovered and numerically checked the following 
interpolation formula [52] 



where the scale 5 g = (1/3) 9 is chosen to be the smallest geometrical feature of the g- 
th generation, and c q and a q are parameters depending on q. This means that for a 
moderately large generation order g, D q! s g , g is almost a linear function oil/g, up to expo- 
nentially small corrections. A linear regression of the values D q ^ g ^ for several generation 
orders allows for a very accurate computation of the limiting values D q of the multifractal 
dimensions. The accuracy of the whole computational method was verified by calculat- 
ing the information dimension D\ which is equal to 1 for any planar simply connected 
set (a result known as Makarov's theorem) [57]. For the quadratic von Koch curve, we 
obtained D\ = 1.0000 ± 0.0001, such a precision being classified as very high for this kind 
of numerical computation. 

2.2.3 Harmonic measure in 3D 

The second important result of [52] is the numerical computation of the multifractal di- 
mensions for the cubic von Koch surface of fractal dimension Dq = In 13/ In 3. For the 
first time, the information dimension D\ of a fractal surface was accurately calculated. 
Contrarily to an intuitive (but false) extension of the Makarov's theorem in three dimen- 
sions, we obtained D\ = 2.007 ± 0.002, providing a concrete numerical counter-example 
to such an extension. It is worth noting that the first mathematical counter-example, for 
which Di was rigorously shown to be strictly greater than 2, was given in [58], but without 
providing the value of D\. Given that the obtained value 2.007 is very close to 2 for the 
cubic von Koch surface, one may wonder whether may exist or not irregular surfaces for 
which the excess D\ — 2 is relatively large. Such surfaces would be particularly promising 
for designing efficient exchangers for diffusion-limited transport. 
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Figure 7: Fourth generation of the self-similar triangular von Koch curve with Hausdorff 
dimension Dq ~ 1.699. Depending on the side exposed to diffusing particles, this curve 
allows one to model either branched pore networks (source at the "top") or fjord-like 
rough bays (source at the "bottom"). 

2.2.4 What makes a boundary less accessible? 

Although the scaling properties of the harmonic measure were largely studied, the re- 
lation between the geometry itself and its multifractal dimensions is still obscure. Is a 
more irregular boundary (greater fractal dimension Dq) more screened? How does the 
presence of deep fjords or a pore network modify the surface accessibility? More generally, 
what makes a boundary less accessible? This last question, which is closely related to 
different optimization problems in chemical engineering, has been addressed in [59]. By 
implementing the GAFRW algorithm, we calculated the multifractal dimensions D q for 
two families of self-similar triangular von Koch curves of variable fractal dimension D 
(Fig. 6). Changing the angle a between two intermediate segments of the generator from 
7r to 0, the fractal dimension Dq of these curves can be varied continuously from 1 to 2: 

~ ln(2 + 2sin(a/2))' 

Depending on the side which is exposed to diffusing particles, the shape of these curves 
is geometrically very different so that one can speak about two families. Starting from a 
distant source at the top of Fig. 7, the particles progressively penetrate into smaller and 
smaller pores of the material. Such curves (called "top-seen") could mimic the geometry 
of a branched pore network. In contrast, the diffusing particles started at the bottom 
of Fig. 7 arrive onto a rough surface with a fjord-like pore structure. The curves of this 
family are called "bottom-seen" . 

The first and quite surprising result is that the multifractal dimensions D q (shown 
on Fig. 8) turn out to be almost identical for the "top-seen" and "bottom-seen" curves 
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Figure 8: Positive-order multifractal dimensions D q of the harmonic measure on the "top- 
seen" (solid symbols) and "bottom-seen" (open symbols) triangular von Koch curves of 
different fractal dimension D : q = 2 (stars), q = 3 (squares), q = 4 (diamonds), q = 5 
(up-pointing triangles), q — 10 (down-pointing triangles), and q = oo (circles). Dotted 
lines are a guide to the eye. 

of the same fractal dimension when Dq < 1.3. Being apparently different geometrically, 
these two types of morphologies are essentially indistinguishable for diffusing particles, 
i.e., when they are "seen" by the harmonic measure. At the same time, the spatial 
distributions of hitting probabilities Pk,s, g over these two boundaries are very different. 
From a certain value of Dq, the multifractal dimensions for the "top-seen" and "bottom- 
seen" curves split. In the former case, these dimensions approach 1, while in the latter 
case, they converge to smaller values. This behavior can be qualitatively explained by 
geometrical arguments (see [59] for details). 

From this study, we draw attention to two facts. First, the fractal dimension Dq, as 
a classical "measure" of the geometrical complexity, is not determinant for the diffusion 
screening: the harmonic measure on the "top-seen" curves with Dq close to 2 exhibits 
almost the same scaling behavior as the harmonic measure on a smooth boundary since 
all the D q are close to 1. This result provides a new insight into practical applications, e.g., 
in chemical engineering. Indeed, it shows that one does not need to explore very irregular 
shapes to diminish the multifractal dimensions. Second, the multifractal dimensions of 
the harmonic measure on the "top-seen" and "bottom-seen" curves are almost identical 
for moderate values of Dq. This means that the harmonic measure is not sensitive to 
distinguish these boundaries of quite different geometry. The hierarchical self-similar 
structure of these curves is thus more determinant than their geometrical details. 

2.3 Spread harmonic measures 

The harmonic measure characterizes only the accessibility of a surface: how the particles 
reach the surface for the first time. When the surface is fully absorbing (infinite reaction 
or transfer rate W), the motion stops immediately after the first hit. When W is finite, 
the particle can be reflected from the surface to the bulk to resume its diffusion. In this 
case, the first arrival point and the final absorption (or reaction) point do not necessar- 
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ily coincide. It is the distribution of the absorption points that is really important for 
transport phenomena. 

2.3.1 Effect of spreading and exploration length 

At microscopic level, a finite transfer rate can be incorporated by introducing a sticking 
probability a: once a particle reaches the boundary, it is either transferred (absorbed, 
relaxed, reacted, etc.) with probability a, or reflected into the bulk with probability 1 — a 
to resume its diffusive motion. For random walks on a lattice with mesh a, the sticking 
probability a was related to the exploration length D/W as a = a/ (a + D/W) [60]. The 
smaller the sticking probability, the larger the number of reflections that corresponds to 
"easy access, difficult transfer" situation (D/W is large). In this case, the particle has to 
explore a certain region of the surface around the first hitting point before being finally 
absorbed (or transferred). While the probability distribution of the arrival points is in 
general sharp and uneven due to diffusion screening, the distribution of the absorption 
points, at which the particles are transferred across the boundary, is getting smoother 
and spread by multiple reflections [61]. Figuratively speaking, reflections "fight" against 
diffusion screening and reduce its effect. In the ultimate case of fully reflecting boundary 
(W = and a = 0), each particle would explore the whole confining domain in an even 
manner, without diffusion screening at all. 

The exploration length D/W controls the spreading effect: larger the D/W, further 
the particles spread around the first hitting point until their absorption. For a flat bound- 
ary (straight line in 2D or plane in 3D), we showed that approximately half of particles are 
absorbed by the disk of radius D/W around the first hitting point [62]. The exploration 
length D/W is thus a physical yardstick for quantifying the geometry. 

2.3.2 Scaling properties of the spread harmonic measures 

For a smooth boundary, the "spread harmonic measure" is defined for any (Borel) subset 
of the boundary as the probability for absorption of the partially reflected Brownian 
motion on this subset. This is a family of measures which are naturally parameterized 
by the length D/W, ranging from the harmonic measure at D/W = to the Lebesgue 
measure at D/W = oo. Then we studied numerically the scaling properties of the spread 
harmonic measures on finite generations of the quadratic von Koch curve [63]. 

The GAFRW algorithm was well suited for simulating Brownian trajectories in the 
bulk. The new implemented feature was partial reflections from the boundary. When 
the particle arrives on the boundary, a random number is generated to decide whether 
the particle is absorbed or not. If not, the particle is reflected (jumped at some small 
fixed distance a from the boundary) to resume its diffusive motion. As the exploration 
length is typically much larger than a, the sticking probability is very small, leading to 
a large number of reflections. The efficiency and rapidity of the GAFRW algorithm were 
therefore crucial for this study. 

As in Sect. 2.2, the spread harmonic measures can be characterized by a set of mul- 
tifractal exponents showing variations of the moments in Eq. (13) with the scale 5. The 
qualitative arguments for determining their scaling properties relied on comparison be- 
tween several length scales: the smallest and largest geometrical lengths £ and L, the 
scale S, and the exploration length D/W. Knowing for a straight line that the exploration 
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Figure 9: Schematic representation of scale zones of the spread harmonic measure. When 
5 is of order of the smallest geometrical detail £, or smaller (on the left), one is looking at 
the harmonic measure of a linear segment. When S is of order of the total size (diameter) L 
of the curve, or larger (on the right), the approximation of the spread harmonic measures 
is too coarse. The most interesting region £ <C 5 <C L includes the transition between the 
Hausdorff and harmonic measures. 

length characterizes the region around the first hitting point where a half of particles are 
absorbed, we conjectured that, for irregular curves, D/W represents a perimeter of this 
region. For a fractal curve, the diameter of this region, £diam, is related to its perimeter 
by a scaling relation, £diam ~ l[(D /W) / £} l l D ° , D being the fractal dimension and I the 
smallest lengthscale of the boundary. The cases £ diam < i and £diam > L trivially cor- 
respond to (almost) harmonic measure and (almost) Hausdorff measure (the Hausdorff 
measure is an extension of the classical Lebesgue measure to fractal objects, see [63] for 
details). The intermediate case with t <C £diam <C L is the most interesting. Depend- 
ing on the ratio 5/^diam, the spread harmonic measures scale differently as schematically 
illustrated in Fig. 9. Figure 10 shows a continuous transition from the multifractal dimen- 
sions of the harmonic measure, D q (studied in Sect. 2.2), to the multifractal dimension 
of the Hausdorff measure, Dq. In addition, the developed concepts brought us to an al- 
ternative derivation of the anomalous constant phase angle (CPA) frequency behavior in 
electrochemistry [42, 64-67]. This new insight allowed us to explain some disagreements 
concerning the CPA exponent (see [63] for details). 

From a mathematical perspective, a continuous transition between the harmonic and 
Hausdorff measures is a new interesting phenomenon. Although the presented study 
remained qualitative and relied on numerical analysis of prefractal curves, the observed 
results suggested a natural way for extending the spread harmonic measures to really frac- 
tal boundaries. The fact that a physical parameter D/W "tunes" the geometrical scale, 
at which the spread harmonic measures are looked at, may have potential applications. 
For instance, the diffusion screening, which is a fundamental obstacle in designing highly 
efficient exchangers, is progressively reduced by increasing the exploration length D/W . 
Moreover, the "functioning" of a system depends on the value of D/W. This means 
that such exchangers can be selective with respect to different species. For instance, this 
principle can potentially be used for designing a filter that would efficiently capture the 
species with a specific diffusion coefficient (i.e., specific size). 
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Figure 10: The scale-dependent multifractal dimension (with q — 2) of the spread har- 
monic measures on generations g — 5, 6, 7, 8 of the quadratic von Koch curve with different 
values of the exploration length D/W. All these dependencies, considered as functions of 
the scaling parameter £ = ln(5/^ diam ), fall onto the same master curve r 2 (£). For small 
£ or 5 (on the left), one retrieves the multifractal dimension of the Hausdorff measure, 
which is equal to the fractal dimension of the curve, D Q = In 5/ In 3 ~ 1.465. For large £ 
or 5 (on the right), r 2 (£) approaches the correlation dimension of the harmonic measure, 
D 2 ~ 0.8925 that we found in [52]. 

2.4 First passage statistics on fractal boundaries 

Reflected Brownian motion is an intermittent process when diffusion steps in the bulk 
are altered by encounters with the surface. In the previous subsection, the active sites, 
which could absorb or transfer the diffusing particles, were supposed to be distributed 
uniformly over the boundary (spatially uniform sticking probability). In this case, the 
overall reaction process was shown to be realized in a region around the first hitting 
point and of the size of the exploration length. When, on the opposite, the active sites 
for reaction or transfer are diluted on the boundary, the surface is not homogeneously 
reactive any more. A more detailed study of individual diffusion steps in the bulk, known 
as Brownian excursions or bridges, is then needed [68]. 

2.4.1 First passage statistics 

For a given boundary, we are interested in the two probability densities (Fig. 11): 

• ip(t) that a particle, started from a close vicinity of the interface at t — 0, returns 
to this interface, for the first time, at time between t and t + dt; 

• and 9(r) that the end-to-end Euclidean distance r between the starting and hitting 
points of the corresponding Brownian excursion (or bridge) lies between r and r + dr. 

For flat surfaces, both densities are known explicitly. For instance, one has for a straight 
line in the plane [1]: 



a 





26 



Figure 11: First passage statistics for a flat surface. Brownian motion started from a point 
S above the boundary (at distance a) hits this boundary at some (random) point H. The 
probability densities ip(t) and 9{r) describe the first hitting time t and the Euclidean 
end-to-end displacement r, respectively 

where a is the height of the starting point above the flat boundary (horizontal axis). 2 
The asymptotic behavior of these densities at long time and large distance is described 
by power laws: 

i)(t)(xt~ a (t-xjo), 0(r)xr~ (r ^ oo), (14) 

with the exponents a = 3/2 and (3 = 2. The first moment (average duration or displace- 
ment of an excursion) and the second moment of the former probability density functions 
are ill-defined and diverge mathematically speaking. But, in most practical situations, 
irregular interfaces are encountered and the behavior associated with flat interfaces is 
possibly misleading. It is thus important to know how an irregular geometry of the in- 
terface may affect the asymptotic behavior (14). In particular, what are the values of the 
exponents a and (3 for a fractal boundary? These issues had been addressed in [68]. 

2.4.2 Scaling relations 

We first derive two general expressions describing the first-passage statistics of Brownian 
excursions (or bridges). A simple qualitative argument is the following. The total number 
Q(t) of particles which have diffused out from the source is approximately proportional 
to the Minkowski content of the surface. The size of this content evolves in time as t 1//2 , 
yielding Q(t) oc [t 1 / 2 ] d_D ° for a boundary of fractal dimension Do, d being the dimension 
of the embedding space [69]. On the other hand, the total number Q(t) is constituted of 
all the particles survived up to time t, Q{t) oc J Q S(t')dt', while the survival probability 
S(t) is related to the first hitting time density ij)(t) as S(t) = 1 — f*i(j(t')dt'. Bringing 

2 It is important to stress that both densities depend on the parameter a. In particular, if a goes to 
0, the densities ?/>(£) and 9(r) approach Dirac distribution S(t) and S(r), respectively. This reflects the 
well known fact that Brownian motion started from a smooth boundary crosses this boundary an infinite 
number of times during an infinitesimal time interval. In other words, Brownian motion immediately 
returns to the boundary. At first thought, this mathematical fact may sound as a puzzling paradox: 
Brownian motion simply cannot escape from the boundary. 

From a physical of view, one should not forget that Brownian motion is a mathematical process which 
leads to the same macroscopic description as real microscopic dynamics (Sect. 1.1). In particular, the 
parameter a cannot be smaller than an interaction range between the surface and the diffusing particle 
(in the order of at least a few angstroms). What is important here is that the parameter a, whatever its 
physically limited value, is much smaller than the transport scales (variable r). In what follows, a will 
be considered as a given small parameter. 



27 



these together, one gets ip(t) oc —d 2 Q(t)/dt and 

Dp-d + A 

a = . (15) 

The displacement statistics 9(r) and the time statistics ip(t) are formally related according 
to 



9(r) = J dt i/>(t) 5^/< r 2 (t) >-r^j, 



where < r 2 (t) > is the mean square displacement at time t. Assuming that the mean 
square displacement of the Brownian motion evolves as t in the bulk phase, a change of 
variable t for the delta distribution in the above equation gives, for a fractal boundary, 
6{r) oc l/r 2 "- 1 , and 

(3 = 2a -I p = D -d + 3. (16) 

The above qualitative arguments give an idea why and how the exponents a and (5 are 
related to the fractal dimension D , without pretending for a mathematical rigor. A more 
rigorous derivation for the exponent (3 in 2D is sketched in [68]. 



2.4.3 Numerical verification 

The scaling relations (15, 16) were checked by extensive numerical simulations for a class of 
self-similar and self-amne interfaces in two and three dimensions. In particular, diffusion 
near various triangular von Koch curves was simulated by the GAFRW algorithm. The 
starting points were chosen uniformly over a finite generation of the curve, within a small 
distance a. The simulation of a Brownian trajectory was terminated when the particle 
approached the boundary closer than a chosen threshold. For each particle, the duration 
and the end-to-end Euclidean displacement were recorded, providing large statistics for 
ip{t) and 9(r) (the number of simulated trajectories was 10 10 ). The data were fitted 
by power laws in order to determine the exponents a and (3. The fact that the fractal 
dimension of the triangular von Koch curves can be continuously varied from 1 to 2 
allowed us for a careful check of the scaling relations (15, 16). A good agreement between 
theoretical predictions and numerical results is shown in Fig. 12. 

A simple dependence of a and (3 on the fractal dimension is a nontrivial result because 
diffusion screening and the harmonic measure properties might influence these exponents. 
The simulations confirmed that the Brownian dynamics in the bulk is not biased by a 
geometrical confinement induced by the boundaries. The exponent f3 is found to be strictly 
larger than 2. Consequently, the mean distance for the first-passage encounter is now finite 
in opposition with the case of a flat surface. The fact that Brownian bridges are sensitive 
to surface geometrical crossovers at long time should provide a way to probe colloidal 
shapes [70]. The studied first-passage process plays a central role in thermodynamics 
of rough colloidal surfaces [71], or in the evaluation of the mean first exit time from a 
bounded domain [72]. It is also important in nuclear magnetic relaxation in complex 
fluids and porous media [73]. 



2.5 Passivation processes 

In the previous subsections, the surface activity remained constant in time. In practice, 
the boundary itself can be altered in the course of transport process. The passivation of 
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Figure 12: Variation of the exponent a (left) and f3 (right) with the boundary fractal 
dimension Dq\ self-similar curves in two dimensions (squares); self-similar surfaces in 
three dimensions (circles); self-affine surfaces in three dimensions (triangles). The solid 
line follows Eqs. (15) and (16), respectively. 

the surfaces working under diffusion-limited conditions is a general phenomenon which 
appears in many natural or industrial systems ranging from catalysis [74] to heat transfer 
[75], electrochemistry [76], and physiology [38]. In such situations, it is not only diffusion 
screening which leads to a strongly inhomogeneous active surface, but part of the surface 
activity may be progressively inhibited by phenomena bearing different names, depending 
on the application field: passivation, fouling, poisoning or restricted absorption. 

For instance, the transfer of nutrients from the digestive system to the blood in humans 
is mostly realized in the small intestine in which the major transport mechanisms are 
passive diffusion and absorption [77]. The anatomy of the small intestine exhibits a 
fractal-like geometry, with finger-like structures at many different scales of magnification: 
flexures, plicae, villi, microvilli [78]. In this type of geometry, the most exposed parts of 
the intestinal membrane are easily accessed by diffusion and thus are the first to be altered 
by any inflammatory disorder or any chemical species that would diffuse in the digestive 
system. As a matter of fact, a wide range of gastrointestinal disorders are associated with 
abnormal intestinal permeability [79]. 

Another important example of passivation is a parallel or serial "fouling" in hetero- 
geneous catalysis [80]. During catalysis, this phenomenon consists in a parasitic reaction 
that passivates the catalyst in the regions which are active, and eventually eliminates 
the entire activity of these regions. Diffusion screening in irregular catalytic grains im- 
plies that active regions represent only a fraction of the total catalytic surface. The time 
evolution of the overall catalytic efficiency will then depend in a complex way on the 
accessibility of the more remote regions of the interface. 

Finally, an enhanced efficiency of heat exchangers is often achieved by building inter- 
faces of very large surface [81, 82]. But the functioning of these interfaces can be sub- 
stantially altered by a fouling process, namely the scale deposition, in which crystalline 
deposits of low thermal conductivity locally reduce the heat transfer [83] . 
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Figure 13: Schematic view of the successive active then passivated regions. The diffusing 
particles are coming from the bottom. The first active region is in red and has a length 
L act . After the first step of the passivation process, the boundary condition on this region 
is set to Neumann, allowing a new region (in green) to become active. This region will 
be in turn passivated and so on, until the whole developed surface is passivated. In 2D, 
the active length L act remains almost constant from one iteration to the next, until the 
entire surface is passivated. 

2.5.1 Iterative passivation process 

In this Section, we investigate the process of progressive passivation of irregular surfaces 
accessed by diffusion [84]. Passivation, disease, aging or fouling will likely start by damag- 
ing the most accessible part of the active interface, so the question then naturally arises: 
what happens after the initially active regions have been passivated? After deactivation, 
most of the diffusing particles hit a now passivated zone and are reflected, resuming their 
diffusion in the bulk to eventually react on an alive but deeper region of the catalytic 
surface (Fig. 13). Consequently, regions that were initially poorly active become fully 
active until they are in turn passivated. This passivation process goes on and on until, 
finally, the whole catalytic interface is deactivated and the catalytic process stops. The 
question of interest here is: how will the size of the region, where most of the particles 
finally react or are absorbed, evolve as the passivation process gradually deactivates the 
initial alive interface. In other words, how the activity is gradually transfered from the 
most accessible to less accessible regions. 

In mathematical terms, the passivation process can be described as follows: at the 
beginning of the process, the alive sites are supposed to be uniformly distributed over 
the whole irregular surface. We suppose that the sticking probability a is equal to 1 
(or W = oo) that corresponds to a homogeneous Dirichlet boundary condition on the 
concentration of reactant molecules. On such an interface, the activity, although existing 
in principle everywhere, is distributed in a very uneven manner due to diffusion screening. 
One may then define an "active zone" as the smallest part of the interface carrying a given 
(large) fraction p of the activity, e.g., 80%. The passivation process is then discretized 
and divided into the following steps [85]: 

• At first, the entire interface is alive. The distribution of arrival probabilities at the 
interface is calculated by the GAFRW algorithm on the quadratic and cubic von 
Koch boundaries as in Sect. 2.2. 
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Figure 14: Successive active regions during the passivation process of a 5th generation of 
the cubic von Koch surface. Diffusing particles are coming from the bottom and reach 
first mainly the blue region of the interface. Each color in the simulation represents a set 
of 4 successive passivated regions (dark blue = regions 1 to 4, light blue = 5 to 8,...). One 
can see that the size of the active region decreases during the passivation process. At the 
end, only the dark red regions on the tip are active. 

• The active region of the interface, which is only a fraction of the alive interface, is 
determined. The activity is proportional here to the harmonic measure density on 
the interface. 

• This first active region is passivated. In mathematical terms, it would correspond to 
a particle concentration obeying Neumann boundary condition. Physically, it means 
that when a particle now hits a site belonging to this passivated region, it is reflected 
back and resumes its bulk diffusion until it reaches the regions that are still alive. 
In other words, this passivation process locally transforms a Dirichlet boundary 
condition into a Neumann boundary condition. The remaining alive interface thus 
consists in the former alive interface minus the newly passivated region. 

• The new distribution of arrival probabilities is now computed using the new bound- 
ary condition (reflecting sites in the former active region). 

• A new active region is determined, which is a subset of the remaining alive interface 
(the non passivated boundary). This new active region will be in turn passivated, 
and so on. 

This passivation process generates at each iteration a new active region so that the whole 
interface can be decomposed as a sum of successive active regions. 

2.5.2 Modification of the GAFRW algorithm 

Although the GAFRW algorithm has proved to be very efficient, its direct application 
to a partially passivated boundary is still not sufficient to solve the problem here. The 
problem is the following: after a few iterations of the passivation process, the remaining 
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alive regions are highly screened. So a reactant particle launched from a distant source 
has to follow a very long stochastic trajectory with a large number of reflections on 
the passivated sites of the boundary before reaching any potentially active region. The 
extremely large computational time required for further passivation iterations makes then 
difficult or even impossible to study the whole passivation process. This difficulty was 
overcome by using the distribution of the activity at step n as the initial source distribution 
for the next step (n + 1). As a matter of fact, the distribution of the activity at step 
(n + 1) is composed of two types of particles: (i) particles arriving directly on the non 
passivated part of the boundary (corresponding to its contribution to the initial harmonic 
measure) and (ii) particles arriving after being reflected by the already passivated part of 
the boundary up to step n. The distant source of diffusing particles is thus replaced by a 
fictitious source on the boundary itself. Since the last passivated regions are close to the 
remaining alive zones, using them as sources considerably enhances the efficiency of the 
computation (see [84] for details). 

This modification of the GAFRW algorithm has permitted us to study the passivation 
of the quadratic von Koch curve (2D) up to the 7th generation and of the cubic von Koch 
surface (3D) up to the 5th generation. One can take note visually of the complexity of 
the 3D surface and the evolution of the active zone in Fig. 14. 

2.5.3 Comparison between 2D and 3D cases 

In 2D, the active region has approximately a constant length L act at each step of the 
process, as schematically indicated in Fig. 13. As one can see in Fig. 15 (left), this 
result is confirmed by numerical simulations for both the 4th and the 7th generations. 
The observed oscillations can be attributed to the discrete scaling of the deterministic 
interface. Moreover, it is known from previous studies that the size of the active region 
is of the order of the width of the interface [57, 85]. This implies that the number of 
passivation steps before the whole interface is deactivated should be of the order of the 
ratio between the total developed perimeter and the size of the interface. Quantitatively, 
the perimeter of the 4th generation (resp. 7th generation) of the quadratic Koch curve is 
(5/3) 4 m 7.7 times (resp. (5/3) 7 w 35.7 times) larger than the size of the system. Hence, 
one can see in Fig. 15 that the length of the active region becomes smaller than 50% of 
the width of the cell respectively after 8 and 37 passivation iterations. After that, in both 
cases the size of the active region falls rapidly. 

In 3D, the striking result is that, unlike in 2D, the surface area of the active region 
S'act is not constant during the passivation process, but gradually decreases (as shown in 
Fig. 15 (right) for p = 80%). Even more, for the 5th generation, the developed surface in 
our 3D simulation contains (13/9) 5 ~ 6.3 times the projected surface, but 25 passivation 
iterations are necessary to completely passivate the surface. We can observe here a net 
discrepancy between the 2D and the 3D cases: due to the properties of Brownian motion, 
the passivation process is much steadier in 2D than in 3D, and it stops much more abruptly. 

This observation is significant for the case of catalysis, where the need of a large surface 
within a finite volume implies the use of very irregular surfaces. As we have shown, the 
deep parts of the surface become more and more difficult to reach as the passivation 
gradually progresses, and the yield of reaction may decrease, despite the fact that active 
regions still exist. In other words, even when catalyst grains seem to be exhausted, they 
may still contain a large amount of alive surface. In this case, a large amount of catalyst 
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Figure 15: Comparison between passivation processes in 2D (left) and 3D (right). Left. 
The size L act of the active region is plotted in units of the system width (L 4 and L 7 for 
generations 4 and 7 of the quadratic von Koch curve) at each step of the passivation 
process. Blue circles and red stars represent generations 4 and 7 in 2D (resp., 3 and 5 
in 3D). One can see that the size L act remains almost constant throughout the process, 
and is of the order of the structure width, until the number of iterations reaches the 
ratio between the perimeter of the interface and its width, L per /L. These ratios are 
represented by the vertical dashed lines. After this threshold, the passivation process 
rapidly terminates. Right. The surface area 5 act of the active region is plotted in units 
of the projected surface (S3 and £5 for generations 3 and 5 of the cubic von Koch surface) 
at each step of the passivation process. In contrast to the left plot, the active region at 
each step regularly decreases. It takes more than 25 steps to passivate the total surface at 
5th generation, whereas the total developed surface represents only 6 times the projected 
surface. The vertical lines correspond to the ratio of the total developed surface (5dev,3 
and Sdev,5) on the projected surface (S3 and S5). Note the progressive slow decrease which 
is very different from the 2D case. 
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Figure 16: Five different geometrical structures with the same surface-to-volume ratio: 
(a) cast of a human acinus, (b) branched Kitaoka labyrinth, (c,d) two long channels 
"packed" in the cube, and (e) a disordered porous medium created by random "digging" 
in the cube. The first four domains satisfy the connectivity condition (accessibility from 
the "entry", indicated by an arrow), while the last one does not. The "solid" channels 
represent the volume of the confining medium where a gas diffuses. 

is wasted simply because the alive surface is not easily accessed by 3D diffusion anymore. 
From this point of view, an engineered (2+l)D geometry (with translational invariance) 
would be preferable. 

2.6 Diffusion-weighted imaging of the lungs 

The above studies were focused on various properties of the diffusive transport from a 
distant source towards an irregular boundary. In other words, we mainly considered what 
is the role of a geometrical irregularity of the interface. For this purpose, finite genera- 
tions of self-similar von Koch boundaries were particularly suited. In many situations, 
however, the internal morphological, topological or geometrical structure of the confin- 
ing domain appears to be even more significant. The examples are porous matrices of 
sedimentary rocks and sandstones, interstitial space of human skin, pulmonary acini, to 
name a few. Although the irregularity of the interface is still relevant, the structure of 
the whole diffusion-confining domain (e.g., pore size distribution, pore connectivity, etc.) 
may dominate the overall transport characteristics. In spite of the active research in this 
field, many questions still remain unanswered. This is probably because of a lack of re- 
liable geometrical models that would be at the same time simple enough for performing 
numerical simulations and representative of porous media found in nature and industry. 
As an example, we refer to a recent progress in a quantitative understanding of the dif- 
fusive transport in the lungs [25, 86-88], which became possible thanks to the Kitaoka 
algorithm for generating a model geometry of the human pulmonary acinus [89]. 

2.6.1 The human pulmonary acinus and its modeling 

A human pulmonary acinus, the gas exchange unit in the lungs, has a dichotomic branch- 
ing structure with tortuous channels densely filling a given volume (Fig. 16a). A "Ki- 
taoka acinus" is designed as a three-dimensional labyrinth of channels with square profile 
(Fig. 16b). Being well appropriate for numerical simulations, this simplified geometry 
captures the essential features of the pulmonary acinus: a dichotomic tree having non- 
symmetric branches of random lengths and filling a given volume. Moreover, the total 
surface area and the average length of the branches correspond to those of real acini. 
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We shall not discuss former works concerning oxygen diffusion in the lungs which were 
clarified by simulations in the Kitaoka acinus [25, 86-88]. In turn, we present Monte 
Carlo simulations that we developed to investigate restricted diffusion of hyperpolarized 
gases (like helium-3) and the consequent signal attenuation in a diffusion-weighted NMR 
experiment [90]. The aim of this study was to facilitate the development of a reliable MRI 
diagnosis of early stage emphysema which results in partial destruction of the alveolar 
tissue and enlargement of the distal airspaces. The geometry of emphysematous acini was 
modeled by removing randomly a fraction of the internal walls from previously generated 
Kitaoka labyrinth. We demonstrated that diffusion-weighted NMR could be sensitive to 
destruction of the branched structure. In fact, partial removal of the interalveolar tissue 
creates "loops" in the tree-like acinar architecture that enhance diffusive motion and the 
consequent signal attenuation. 



2.6.2 Monte Carlo simulations accounting for magnetic field gradients 

Given the specific shape of the Kitaoka acinus (long branched channels of the same square 
profile), there was no need for using fast random walks. An accurate accounting for a 
linear magnetic field gradient would further complicate an implementation of such an 
algorithm. For this reason, we limited ourselves to implementation of simple off-lattice 
random walks with normal reflections on the boundary of the Kitaoka acinus [90, 91]. 
The given temporal profile fit) of the applied magnetic field gradient was discretized on 
the simulation time interval [0,t] with a time step r = t/n. The starting point ro is 
chosen randomly with a uniform distribution inside the Kitaoka acinus. For each step 
k, one generates independent Gaussian displacements dr l in the three space directions 
(i = x, y, z) with mean zero and dispersion \ / 2Dt, in order to pass from the current 
position r k to a new position r k+1 . If the linear segment between the current position 
and the new position intersects the boundary, a mirror reflection is applied. At each 
step k, the term Tf(kr)r l k is added to the phase counter (p^ for each space direction 
i — x,y, z. The total phase accumulated during the whole trajectory for unit gradient is 
then approximated, for each i, by the sum 

n 

If the gradient is applied along a given direction e = (e x , e y , e z ), the total phase accumu- 
lated in this direction for unit gradient is ip e = e x ip x + e y ip y + e z ip z . Although the gradient 
direction is fixed in the NMR scanner, the real acini in the lungs are oriented randomly. 
This effect is taken into account by looking at the signal that is averaged over all spatial 
directions. In this case, the averaged phase is ip av = y/ip x + + vl- Repeating the 
Monte Carlo simulation iV times, one records these phases in order to approximate their 
probability distribution. Once the simulations are terminated, one can find the signal E 
or the directionally averaged signal S ay (g) as a function of the gradient amplitude g as 

Since the statistical error of Monte Carlo technique is of order of 1/y/N, the choice of 
N = 10 6 leads to reasonably accurate results. 
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Figure 17: Directionally averaged signal S^ig) as a function of the gradient intensity g 
for different destruction factors v starting from (healthy acinus) and ranging up to 1 
to model progressively damaged emphysematous acini. The signal is more attenuated in 
more damaged structures because diffusion is faster with a lower restriction. The ratio 
between the signals from a healthy acinus iy = 0) and that from an early emphysematous 
acinus [y = 1/6) is higher for larger gradient intensities. 

2.6.3 Diffusion in healthy and emphysematous acini 

Figure 17 illustrates the dependence of the directionally averaged signal S av (g) on the 
gradient intensity g for healthy and emphysematous acini. For weak gradient intensities, 
one recovers the classical (Gaussian) g 2 behavior of the logarithm of the signal, with a 
single apparent diffusion coefficient (ADC). Since partial destruction of the alveolar tissue 
by emphysema creates loops in the branched structure, diffusion becomes faster, and the 
signal is then more attenuated (larger ADC). The high sensitivity of diffusion- weighted 
NMR measurements to this effect can be potentially employed to diagnose emphysema at 
early stages. 

The geometrical confinement and branching structure of the acinus lead to deviations 
from the g 2 behavior at higher gradients. In particular, one can see on Fig. 17 a transition 
to a stretched-exponential behavior known as localization regime [92, 93]. This observa- 
tion indicates that the notion and use of ADC should be substantially revised, especially 
because experimental measurements at higher gradients appear to be more sensitive to 
the acinar structure. This also explains several confusions and possible ambiguities in 
determination of ADC in medical literature. 

The numerical simulations were also performed on a tenfold scale Kitaoka acinus 
with a broad range of diffusion coefficients. The obtained results were confronted to 
experimental measurements in a tenfold phantom made from the standard epoxy resin by 
stereolithography [94, 95]. 

2.6.4 What is the role of a complex internal architecture? 

In a separate work [96], we addressed a more general question: What is the role of a 
complex internal architecture (e.g., branching of the pulmonary acinus or pore network 
in rocks) for NMR measurements? To answer this question, we have performed Monte 
Carlo simulations of restricted diffusion in three groups of three-dimensional structures 
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Figure 18: Directionally averaged signal S av (g) as a function of the gradient intensity g for 
different porous structures on Fig. 16 (the domains called "disordered-2" and "disordered- 
3" are not shown). 



with the same surface-to-volume ratio. A basic domain is a cube of size L divided into 
6x6x6 = 2 16 small cubic cells. The first group is a set of random dichotomic labyrinths 
generated inside the cube by the Kitaoka algorithm (Fig. 16b). The second group con- 
sists of two realizations of a long channel filling the same cube (Fig. 16c, d). In the third 
group, disordered porous media are generated by connecting a number of randomly cho- 
sen adjacent cells (Fig. 16e). The directionally averaged signals from these structures are 
computed by the above Monte Carlo simulations. Figure 18 shows that the signal atten- 
uation Sa^g) for different branched structures (the first and second groups) are almost 
identical, while the signal for disordered media (the third group) is significantly higher. 
This is related to the fact that a disordered medium consists of a number of small dis- 
connected patterns where the signal is less attenuated. The important conclusion is that 
the internal structure of porous media significantly influences the restricted diffusion and 
NMR measurements and should thus be taken into account for practical applications. 



Summary 

In this Section, we showed that Monte Carlo simulations were flexible and efficient tools for 
studying restricted diffusion in complex geometries, e.g., von Koch fractals. A geometry- 
adapted fast random walk algorithm allowed us to solve different problems, including 
multifractal properties of the harmonic measure in 2D and 3D (Sect. 2.2), scaling prop- 
erties of the spread harmonic measures and the role of the exploration length (Sect. 2.3), 
first passage statistics (Sect. 2.4), passivation processes (Sect. 2.5), etc. All these prob- 
lems concern restricted diffusion with different behaviors on the boundary. The GAFRW 
algorithm simulated Brownian trajectories in the bulk, while the boundary behavior was 
specifically implemented for each problem. 
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3 Spectral insight: Laplacian eigenfunctions 



This section is devoted to the spectral description of restricted diffusion when the Laplace 
operator eigenfunctions allow one to structure the whole information about diffusion in 
a particularly useful form. Most diffusion characteristics (e.g., propagator, total flux, 
residence time, etc.) can be explicitly written in terms of the eigenfunctions. Specific 
features of these characteristics originate somehow from eigenfunctions. We aim to reveal 
their relation and understand various features of restricted diffusion in complex geometries 
via the properties of the underlying Laplace operator eigenfunctions. While this spectral 
description appears to be natural, it is not always considered as the most efficient way. 
In fact, except few simple domains, for which the Laplace operator eigenfunctions are 
known explicitly, a numerical calculation is required. This is a difficult time-consuming 
computational problem, especially for porous, multiscale, or irregular geometries. But, 
as far as the eigenfunctions are found for a given domain, the whole range of diffusion 
characteristics becomes directly accessible at once. 

We used the spectral description as a unifying mathematical language for presenting 
the main achievements on restricted diffusion in NMR [14, 118, 119]. The diversity and 
complexity of diffusive NMR phenomena, observed in experiments, were shown to result 
from the specific properties of the reflected Brownian motion and the underlying Laplace 
operator eigenfunctions. Many classical results were retrieved, extended and critically 
discussed. In what follows, we present our main results basing on this spectral point of 
view. 

3.1 Matrix formalism 

Several matrix formalisms were developed for numerical analysis of the macroscopic signal 
formed by the nuclei diffusing in magnetic fields [27, 120-123]. We have reformulated and 
extended these formalisms to describe restricted diffusion in any geometrical confinement 
and arbitrary magnetic field [14, 124]. The "derivation" of the compact matrix form for 
the signal (shown in the next subsection) simply relies on a representation of the Bloch- 
Torrey equation (4) in the Laplace operator eigenbasis. The use of this natural basis 
allows one to get a structured representation of diffusion via two governing matrices A 
and B. A truncation of these matrices for further numerical analysis is well controlled, 
leading to negligible computational errors. The use of the Laplace operator eigenfunctions 
is crucial here. 

3.1.1 Perturbative derivation 

When the applied magnetic field is independent of time, B(r,t) = /3B(r), Eq. (4) reads 
as 



where B(r) is the normalized (dimensionless) spatial profile of an inhomogeneous magnetic 
field, and (3 its intensity. We are looking for a solution of this equation in the basis of the 
Laplace operator eigenfunctions u m (v) (of eigenvalues A m ) satisfying Eqs. (9a, 9b): 




(17) 




(18) 
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with unknown time-dependent coefficients c m r(t). By construction, c(r,t) satisfies an 
appropriate boundary condition that was imposed for the eigenfunctions. Substitution 
of this expansion in Eq. (17), multiplication by -u^(r), and integration over the confining 
domain fl yield a set of ordinary differential equations on the coefficients c m (t) 



dc m (t) 
dt 



+ J2\ DAm > m ' +n(3Bm,m')cm'(t) = 0, (19) 



where the infinite-dimensional matrices A and B represent the Laplace operator and the 
"perturbing" magnetic field B(r) in the eigenbasis of the Laplace operator: 



J dr <(r) B(r) u m ,(r), (20) 



A m ,m' — $m,m'^m- (21) 

Thinking of c m (t) as components of an infinite-dimensional vector C(t) leads to a matrix 
first-order differential equation 

j t + DA + i-fpB^j C(t) = 0. (22) 

As for a scalar equation, its solution is the (matrix) exponential: 

WC(t) = Ue- iDA+ ^ B)t . (23) 

The supplementary factor \/V is put explicitly to get rid of the dimensional unit, meter" '/ 2 , 
of the vector C(i) (V being the volume of the domain). Here the matrix exponential 
e -(DA+t-yi3t3)t ac ^ g on i e ft on the vector U representing the initial density p(r) in the 
basis {u m (r)}: 

U m = V 1 ' 2 J dvu* m (v) p(v). (24) 
n 

The macroscopic signal is then obtained according to Eq. (12) by integrating the trans- 
verse magnetization c(r, t) over the whole confining domain Q with a sampling or pickup 
function p(r) of the measuring coil or antenna: 

E = J dr c(r, t) p(r) = ]T c m (t) V^ 2 J dr u m (r) ~p(r) . (25) 



The last sum can be interpreted as a scalar product between the vector C(t) and the 
vector U representing the pickup function p(r) in the basis of the eigenfunctions {u m {r)}: 

U m = V- 1 ' 2 J dvu m {v)~p{v). (26) 
n 

The macroscopic signal at time t can thus be written in a compact matrix form of a scalar 
product: 

e= ( Ue -m+wB)t fj^ (27) 
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matrices being acting on the left. From a quantum-mechanical point of view, the matrix 
e -(DA+i 7 /3B)t can be thought of as an evolution operator acting on the initial state p(r) 
(represented by the vector U), while the resulting density c(r, t) is weighted by the pickup 
or sampling function p(r) (represented by vector U). A rapid increase of the eigenvalues 
\ m with m allows one to truncate the matrices A and B to moderate sizes in order to 
perform the numerical computation. 

Although the compact matrix form (27) is derived for time-independent field B(r), 
the approach is easily applicable to various NMR sequences. For instance, Carr-Purcell- 
Meiboom-Gill (CPMG) sequence consists in repeating the inverting 180° radio-frequency 
pulse n times to acquire successive echoes [125]. In the matrix formalism, the amplitude 
of the fcth echo (k = 1, ...,n) is simply calculated as 



E k = U 



e -{DA+vrpB)t/ (2n) e -{DA-vrpB)t/ (2n) 



k 

u 



We developed a spectral analysis of the multiple echo attenuation for CPMG sequences 
[126]. 

In general, one can approximate any given temporal profile f(t) by a piecewise constant 
function: f(t) = fa on [tk, ifc+i] (k = 0, K). For each interval [tk, tk+i], the Bloch-Torrey 
equation can be solved using the matrix exponential. In order to merge the solutions, 
the ending magnetization of the interval [tk-i,tk\ is taken as the initial condition for 
solving the problem on the next interval [tfc,£fc+i]. The use of the matrix representation 
is particularly efficient to handle such successive computations: 



E=[U 



■ K 

n 

■k=0 



e -(DA+ij/3f k B)(t k+1 -t k ) 



U), (28) 



where to — and t K+ i = t. It is worth stressing that this result is exact, no approximation 
was involved. Moreover, since any function can be approximated by a piecewise-constant 
function, the above relation allows one to approximately compute the signal for arbitrary 
temporal profile f(t). 

Conceptually, Eq. (27) is nothing more than a representation of the Bloch-Torrey 
equation in the Laplace operator eigenbasis. The matrices B and A, depending merely 
on the eigenbasis and the spatial profile B(r), have to be computed only once for a given 
geometry. Once this preliminary step is achieved, the computation of the signal for a given 
set of physical parameters is straightforward, accurate and very rapid. This is the crucial 
advantage with respect to conventional numerical techniques for solving the Bloch-Torrey 
equation (Monte Carlo simulations, finite difference or finite element methods, etc.). 

As we pointed out in Sect. 1.4, the computational efficiency of the matrix formalism 
should not be surprising because finding the eigenbasis is equivalent to solving all the 
diffusive problems at once for a given geometry. Naturally, finding the eigenfunctions 
is a difficult task, especially for complex geometries. But once the eigenbasis is found, 
the remaining computations are easy. We are going to discuss the recent progress that 
we achieved by using the matrix formalism for simple geometries, for which the Laplace 
operator eigenbasis is known explicitly: a slab (an interval), a cylinder (a disk), and a 
sphere, as well as circular and spherical layers. 
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3.1.2 Moments 



The compact matrix form (27) of the signal is particularly suitable for numerical pur- 
poses because the computation of matrix exponentials is rapid and very accurate. At the 
same time, a theoretical analysis with this form faces considerable difficulties because the 
governing matrices B and A do not commute. In this subsection, we follow Ref. [14] and 
briefly discuss another strategy, when the signal is represented as a power series with the 
moments of the normalized total phase <f> = ip/t: 

E = E{e^} = jr ^^E{0 ra }, (29) 

n=0 U ' 

where q = jflt is a dimensionless parameter characterizing the intensity of the applied 
field. Each moment can be explicitly written by using the definition (5) of 0: 

t t t 

W"} = /(*0 / dh /(f 2 )-.. J dt n f(t n ) E{B(X tl )...B(X tn )}. (30) 

tl tn-l 

where the time moments ti, ... , t n were put in ascending order. The multiple correlation 
function E{B(X tl )...B(X tn )} can be calculated using the Markov property of Brownian 
motion: 

E{B(X tl )...B(X tn )} = J dr p(r„) j dn G^n) B(n) ... 

n n 

J dr n G tn _ tn _ 1 (r n _i,r n ) B{r n ) J dr n+1 (7 t _ tn (r n , r n+1 ) p(r n+1 ). 

This lengthy expression has a simple probabilistic interpretation. Starting at time from a 
randomly chosen point r (with a given initial density p(r )), a particle diffuses up to time 
ti to point i - ! where it "experiences" the magnetic field B(ri). The diffusive propagator 
Gr tl (r , i"i) describes this process. From r l5 the particle diffuses to a new point r 2 , and 
so on. At last, the particle arrives at point r n at time t n and experiences the magnetic 
field B(r n ). During the the remaining time t — t n , the particle can diffuse to another r n+1 
where it is "sampled" with a given function p(r n+ i). Using the spectral decomposition of 
the diffusive propagator Gt{r, r'), one can deduce a compact matrix form for the multiple 
correlation function [14] 

E{B(X tl )...B(X tn )} = {Ue- DMl Be- DK ^- tl) B ... # e -^(*n-tn-i) Be -DA(t-t„)^ ^ 

with the same matrices B and A and vectors U and U as earlier. This is simply a 
representation of the above integral form in the Laplace operator eigenbasis. 

In many practical cases, the boundary is purely reflecting (Neumann boundary condi- 
tion), while the initial density p(r) and sampling function p(r) are uniform. It is easy to 
check that the components of the vectors U and U are then zero, except for the ground 
mode: U m = U m = S mt0 . Given that the matrices e~ DAtl and e - DA (*-*«) are diagonal, the 
multiple correlation function reads as 

E{B(X tl )...B(X tn )} = [Be- DA ^-^B ... Be~ DK ^- tn -^B] QQ . 
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For instance, the second moment is 



E{0 2 /2} =< £ B , m e -^(^) Bmfl > 2 , (32) 

m 

where < ... > 2 denotes the time integral in Eq. (30), for instance, 

t t 

<e -DX m (t 2 - tl)> ^ = J dh f ^ e -DA m (t 2 - tl ). 

*i 

In this formal way, the original problem of finding the macroscopic signal of diffusing 
spins is entirely reduced to the analysis of the Laplace operator eigenmodes, and is thus 
solved as a physical problem. In what follows, we show how this mathematical basis can 
be applied for a theoretical analysis in many cases of particular interest. 



3.2 NMR survey of restricted diffusion 

In this section, we present various results on restricted diffusion in magnetic fields that 
we could obtain by using the Laplace operator eigenfunctions [14]. 



3.2.1 Simple geometries 

Restricted diffusion in a slab (an interval), a cylinder (a disk), and a sphere is a classical 
problem [34, 35]. Since the Laplace operator eigenfunctions in these domains are known 
explicitly, many diffusion characteristics can be calculated analytically. For instance, the 
eigenfunctions and eigenvalues of the Laplace operator in the unit disk are 

u nk (r) = — — J n {a nk r) cosmp, X nk = a nk , 

y/iT J n {a nk ) 

where e n = a/2 — <5„ i0 , the constants (3 nk = a/ \ tk / (^nk — n 2 + h 2 ) guarantee the normal- 
ization (10), J n (z) are the Bessel functions of the first kind, and a nk are all the positive 
roots of the equations 

z J' n {z) + hJ n {z) = 0, (33) 

coming from the Robin boundary condition (9b), with h = W/D (L = 1). The fac- 
torization of the dependences on the radial and angular coordinates r and (p is a direct 
consequence of the rotational invariance. It is worth stressing that here, the single index 
m is replaced by a double index nk (with n — 0, 1, 2, 3... and k — 0, 1, 2, 3...) coming from 
enumeration of the positive roots a nk . We shall use the double index as a convenient 
notation to enumerate the eigenvalues, eigenfunctions and the elements of matrices and 
vectors. 

The explicit form of the eigenfunctions is employed to investigate the signal attenua- 
tion due to restricted diffusion in inhomogeneous magnetic fields (see [14] and references 
therein). This knowledge facilitates the construction of the matrix B for a given profile 
B(r). In general, this would require numerical integration in Eq. (20), but for some choices 
of the function B(r), the integration can be realized analytically. In [14], we provided, 
for the first time, the explicit formulas for the matrix B in two cases of practical interest: 



42 



a linear magnetic field gradient, B(r) = x, corresponding to usual experimental setup 
for NMR measurements, and a parabolic magnetic field, B(r) = |r| 2 . For example, we 
obtained for the unit disk in a linear gradient: 

n a n i a la ^1/2/3 a Xnk ± Xn ' k ' Z 2nn> + 2h ( h ~ _j (*A\ 

O n k,n'k' — On,n'±l{± + Onfl + On' ,0) PnkPn'k' 7T 7 ^ • [o*) 

\^nk — \i'k') 

This matrix turned out to be fully expressed in terms of the eigenvalues X n k- This means 
that the originally complicated problem of computing the macroscopic signal attenuation 
is reduced to finding roots of the Bessel functions in Eq. (33). This result allows one for 
efficient numerical computation of the signal, as well as for a thorough theoretical analysis 
of the moments, as illustrated below. 



3.2.2 Time-dependent diffusion coefficient 

Apart from being a basis for efficient numerical computations, the matrix formalism is a 
powerful analytical tool. As a representative example, we consider a linear magnetic field 
gradient in a given direction e and of intensity g: 

(3 = gL, B(r) = (e • r)/L, 

L being the characteristic size of the domain. For relatively small gradients (small q) , the 
macroscopic signal in Eq. (29) is essentially determined by the second moment E{0 2 /2}: 

E~l- g 2 E{0 2 /2} + 0(q A ) ~ exp[-g 2 E{0 2 /2}] + 0(q 4 ) (35) 

(the first and other odd moments that would determine the phase of the signal are omit- 
ted). This formula known as "Gaussian phase approximation" (GPA) was widely em- 
ployed in many theoretical, numerical and experimental studies of restricted diffusion 
(see [14] and references therein). It is worth noting that this formula becomes exact for 
free (unrestricted) diffusion, for which the second moment was calculated by Stejskal and 
Tanner for arbitrary temporal profile f(t) of the applied gradient [127] 

E{0 2 /2}o = ^<(ti-t 2 )> 2 . (36) 

Deviation of the second moment E{0 2 /2} for restricted diffusion from its value E{0 2 /2} o 
for free diffusion can be used to characterize a confining domain. Actually, their ratio is di- 
rectly related to apparent, effective or time-dependent diffusion coefficient D(t) measured 
in NMR experiments 3 

D(t) = D E{02/2> 
V[t) ^E{0 2 /2V 

According to Eq. (35), this quantity is simply proportional to the logarithm of the macro- 
scopic signal and is thus easily accessible in experiments. A number of studies related 
the behavior of the time-dependent diffusion coefficient and the geometry of a diffusion- 
confining domain [97]. 



3 It is worth noting that this NMR definition of the time-dependent diffusion coefficient is different 
from its classical form, the latter being obviously independent of the applied magnetic field (see [14] for 
further discussion). 
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If there is no surface relaxation, Eq. (32) yields a general representation of the time- 
dependent diffusion coefficient for any confining domain and arbitrary temporal profile 
[14, 128] 

^ = Y, L2X mBl m w f {D\ m tl (37) 

m 

where the function u>/(p) accounts for a given temporal profile fit): 

< e ~p(t2-tl)/t > 

Wf (p) = __2_. (38) 

It is easy to check the following normalization property for linear gradients (taking L = 1 
for simplicity): 



5> 



B 2 

- D 0,m 



so that \ m B 2 m can be considered as the weight of the mth eigenmode to the time- 
dependent diffusion coefficient. In turn, the eigenvalue X m determines a characteristic 
scale at which the contribution appears in Eq. (37) through the function Wf(p). It is 
important to stress that Eq. (37) allows one to distinguish the roles of various "ingredi- 
ents" of the problem. So, the spatial profile of the magnetic field enters only into the 
weights A m £>Q m , while the temporal profile is fully taken into account by the function 
Wf(p). Finally, the geometry determines both the time scales (DA m ) _1 and the weights 
A m £>o m . This connects directly the spectral properties of the Laplace operator to mea- 
surable characteristics of complex media. 

The general representation (37) allowed us to retrieve and extend many classical re- 
sults. In the long-time diffusion regime, when the diffusion length \f~Dt is much larger 
than L, the function Wf(p) behaves as 

w f (p)e,p- 2 J^J^L + 0(p- s ) (p»l), (39) 

yielding 

E ~ exp 



7V^ 4 



t 

2/ 



C-i dtf\t) 



o 



D 

for any bounded domain, with the geometry-dependent dimensionless constant 

c-i = £fiU L2A ™) _1 - 

m 

This is an extension of the classical results by Robertson and Neuman [129, 130]. It shows 
a high sensitivity of diffusion-weighted measurements to the size L of a confining domain 
in the long-time regime. 

In the opposite short-time diffusion regime (p 1), only a small fraction of particles 
near the boundary can "feel" their geometrical confinement, while the remaining majority 
of particles diffuse as they were in free (unbounded) space. The fraction of "restricted" 
particles can be estimated as the ratio between the volume of the diffusion layer of width 
\J~T)i near the boundary and the total volume: Sy/Dt/V, S being the total surface area. 
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The deviation of the time-dependent diffusion coefficient from its nominal (free diffusion) 
value D is then quantified by this ratio. A more rigorous analysis based on Eq. (37) gives 
[14] 

m^.m^s <fe-;f >> +oit) , (40, 

D < (*2 - *l) > 2 Vt 



geometry v 

temporal profile 



where C3/2 is a geometry-dependent constant. This is an extension of the classical result 
by Mitra et al. [131-133]. This relation was suggested for experimental measurement of 
the surface-to- volume ratio S/V in porous media. 



3.2.3 Different diffusion regimes 

In many cases of practical interest, the knowledge of the second moment is enough for an 
accurate approximation of the macroscopic signal. When the gradient intensity increases, 
the fourth and higher-order moments become progressively more and more significant (see 
also Sect. 3.2.5). In particular, Stoller and co-workers obtained the asymptotic behavior 
of the macroscopic signal at high gradient intensity for a slab geometry [92]: 

(41) 

where a\ ~ 1.0188 is the absolute value of the first zero of the derivative of the Airy func- 
tion. In this so-called localization regime, the signal is essentially formed by the nuclei 
diffusing (or "localized") near the boundary and thus acquiring less dephasing than the 
nuclei in the bulk. Four years later, this theoretical prediction was experimentally con- 
firmed by Hurlimann et al. [93]. Measuring the signal attenuation from water molecules 
diffusing between two parallel plates, Hurlimann and co-workers observed a spectacular 
deviation from the Gaussian g 2 dependence of InE at gradient intensities higher than 
15 mT/m. 

In our review [14], we summarized different diffusion regimes (Fig. 19). Since diffusion 
and encoding in a spin-echo experiment can be set independently, we plot the signal as 
a function of the two dimensionless parameters p = Dt/L 2 and q = jgLt. From this 
plot, one can determine which kind of restricted diffusion is to be expected, and which 
formula should be applied to fit experimental data. It is worth noting, however, that this 
information is rather qualitative since the particular location of different regions on the 
diagram depends on the diffusion-confining geometry and applied magnetic field. 



E oc exp 



-|(zWt 3 ) 1/3 



3.2.4 Rigorous results for relaxing boundaries 

For three simple domains (interval, disk and sphere), many relevant quantities can be 
found analytically even in the presence of surface relaxation. Rigorous computations 
rely on explicit formulas for the elements of the matrices B and A that were shown to 
be rational functions of the eigenvalues (e.g., Eq. (34)). Using the Laplace transform 
summation technique, we derived exact and explicit representations for the zeroth and 
second moments in the presence of surface relaxation [134]. Within the Gaussian phase 
approximation, these two moments determine the reference and diffusion-weighted signals, 
respectively. 

Without presenting the details, we mention that the computation was essentially based 
on two observations: 
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Figure 19: Different regimes of restricted diffusion in a slab under a linear magnetic field 
gradient. The signal is attenuated by a factor of 2 at each line separating two adjacent 
gray-scale regions (appearing as pale and dark stripes). The first large pale region on 
the left is composed of points (q,p) for which the signal E lies between 1/2 and 1. The 
next dark stripe regroups points (q,p) for which 1/4 < E < 1/2, and so on. The white 
area on the right corresponds to pairs (q,p) for which the signal is below 10~ 3 . Since 
such a small signal is often comparable to noise, this area is referred to as inaccessible 
experimentally. On this pq diagram, the bold line delimits the region on the left in which 
the GPA predictions of the signal attenuation are valid with an accuracy of at least 5%. 



The Laplace transform of the even-order moments can be reduced to a combination 
of multiple sums of the form 

^ 1 1 1 1 1 

J-^ m {?\ ~ (V - A m2 )^ ( S2 _ A m2 )A (A m2 - \ )fa - ( Sn _ \ m y n 

where {sk} is a set of real parameters, and {flk} is a set of positive integer numbers. 

Algebraic transformations and differentiations can reduce any sum of this form to 
the function 



m 



which can be calculated explicitly for the eigenvalues A m of the Laplace operator in 
the interval, disk, or sphere. 

As a consequence, the Laplace transform of the even-order moments have explicit, though 
cumbersome, analytical forms. In particular, these explicit forms allow one to study the 
asymptotic behavior of the moments in different diffusion regimes. In the short-time 
regime, the series expansion in half-integer powers of the diffusion coefficient was gener- 
alized to arbitrary temporal profile of a linear magnetic field gradient. In the long-time 
regime, it was shown how the presence of surface relaxation modified classical Robertson's 
results (for details, see [134]). 
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Figure 20: Numerically computed second moment E{0 2 /2} for a linear gradient with 
bipolar temporal profile as a function of p = Dt/L 2 and its explicit approximation derived 
from the exact solution (43) for a cosine spatial profile. The inset shows relative deviation 
between them. 

3.2.5 Cosine magnetic field in a slab 

Although the above analysis could in principle be carried out for any even-order moment, 
its practical implementation was already tedious even for the second moment. 

But a better knowledge of the high-order moments E{0 n } would help to understand 
the limits of applicability of the Gaussian phase approximation and the transition to- 
wards non-Gaussian regimes. We addressed this problem for restricted diffusion between 
parallel planes in a cosine magnetic field [135]. The specific choice of the spatial profile 
to be proportional to an eigenfunction of the Laplace operator in this confining geom- 
etry considerably simplified the underlying mathematics. In fact, the eigenfunctions of 
the Laplace operator on the unit interval (L = 1) with reflecting endpoints are simply 
u m {x) = e m cos(7rmx), so that the substitution of the spatial profile B(x) = cos(7rx) into 
Eq. (20) yields a particularly simple form 



where e m = a/2 — 5 m ,o are the normalization constants. This subdiagonal structure of 
the matrix B suggests the use of an analogy with reflected random walks to calculate 
the multiple correlation functions E,{B(X tl )...B(X tn )}. For instance, the second moment 
with the cosine magnetic field is simply 



Exact and explicit relations for several higher-order moments were also reported [135]. 

These results helped us to study the general structure and the properties of the high- 
order moments, as well as deviations from the GPA at magnetic fields of moderate or high 
intensity. Although the cosine magnetic field is indeed very specific (e.g., its experimental 



B, 




1 



E{0 2 /2} = ^ < eM-Dn 2 (t 2 - tj] > r 



(43) 
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Figure 21: Circular layers of thickness i = 
L — 1: R = 0.5 (a) and R = 0.9 (b). 
coordinates are orthogonal to each other. 



1 — R with inner radius R and outer radius 
Note that variations in radial and angular 



realization presents a challenge in itself), the deduced properties of the high-order mo- 
ments provide a better understanding of the GPA and its limitations in general. Moreover, 
numerical simulations showed that the explicit relations obtained for the cosine magnetic 
field can be used as good approximations for linear magnetic field gradients (as illustrated 
in Fig. 20 for the second moment). Note also that Zielinski and Sen considered the cosine 
spatial profile to model and study susceptibility-induced local magnetic fields [136]. 

3.2.6 Diffusion in circular and spherical layers 

There still exists a substantial "gap" in understanding restricted diffusion in simple geome- 
tries, on the one hand, and in natural porous media such as sedimentary rocks, cement, or 
biological tissues, on the other hand. From a theoretical point of view, the complexity of 
porous structures is related to multiple length scales, ranging from the size of tiny pores 
(in the order of several microns in some rocks) to the overall size of the sample (e.g., 
few centimeters). Moreover, porous media often present a hierarchical structure of pores 
so that intermediate characteristic length scales emerge. In this case, even the introduc- 
tion of representative dimensionless parameters such as p, q, and h becomes ambiguous 
as being depended on a geometrical length L used. It is thus important to understand 
how classical theories for single-scale shapes (like a slab) can be extended to incorporate 
multiple scales. Since an analytical calculation in natural porous media is practically 
infeasible, some simplified confining domains with two or multiple length scales become 
particularly useful. 

In this light, circular and spherical layers (or shells) appear as appropriate models 
(Fig. 21). Like for a disk and a sphere, the Laplace operator eigenbasis in these domains 
is known explicitly. Using the properties of the underlying Bessel equation, we derived 
explicit analytical formulas for the matrix B representing a linear magnetic-field gradient 
(or parabolic magnetic field as well) in the Laplace operator eigenbasis [128]. From the 
numerical point of view, the problem is completely reduced to finding roots of the equa- 
tions with Bessel functions. Once these roots are found, the governing matrices B and A 
can easily be constructed and used to compute the macroscopic signal for any temporal 
profile of the magnetic field. 

Bearing in mind multiscale structures, we particularly focus on thin layers (Fig. 21b), 
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Figure 22: Normalized time-dependent diffusion coefficient D(t)/D as a function of the 
dimensionless parameter p = Dt/L 2 (L = 1) for circular layers with inner radii R = 0.5 
(left) and R = 0.9 (right). The solid line shows the exact computation via Eq. (44), while 
the dashed line represents an approximate solution containing only the first contributing 
eigenvalue Aio and thus failing for small t. Two vertical dashed lines (at £ 2 and 1) roughly 
split the p axis in three regions: short-time regime (on the left, p -C £ 2 ), intermediate 
region, and long-time regime (on the right, p ^> 1). The vertical dotted lines show the 
scales X^ k (k ranging from to 10, with k = on the right and k = 10 on the left). One 
can observe a plateau due to a large separation between the two contributing eigenvalues 
Aio and An (the thinner the layer, the wider the plateau). Note that 2D diffusion in a 
thick circular layer (a 2D shape) is reduced to essentially one-dimensional motion in a 
thin circular layer, yielding the reduction factor 1/2 (the level of the plateau). 



for which perturbative calculations turned out to be surprisingly accurate (here we assume 
that there is no surface relaxation). For instance, the perturbation series of the eigenvalues 
in powers of the small dimensionless thickness £ are (in 2D): 

A,o = n 2 (l + 1 + \e> + \e - + 0(£ 5 )) , 

^^(i + ^W^ + o^)) <*>o), 

so that the analysis is indeed much simpler for thin layers than for the unit disk (we took 
L = 1 so that the eigenvalues are dimensionless). The smallness of the thickness I with 
respect to the perimeter 2tt results in a large gap between A„o ~ n 2 and A n & ~ 7T 2 k 2 /£ 2 with 
k > 0. The smaller eigenvalues {A„,o} describe large displacements in angular coordinates 
along the boundaries (i.e., along the circumference of the circular layer). In turn, the 
larger eigenvalues {X n k} (with k > 0) describe small displacements in radial coordinate 
between the inner and outer boundaries. This situation is somehow analogous to restricted 
diffusion in a thin rectangle with sides tt and £, for which X™^ = n 2 + 7r 2 k 2 /£ 2 . 

Apart from the above perturbative analysis and derivation of the explicit formulas 
for the matrix B, the main result of [128] concerned the behavior of the time- dependent 
diffusion coefficient in these two-scale domains. Although the summation in Eq. (37) was 
carried out over all eigenmodes (here enumerated by double index nk), the rotation in- 
variance of the circular and spherical layers eliminates the contributions of all eigenmodes 
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except for n = 1: 



D(t) 



D 



X 1Q B 2 10 w f (DX 10 t) + \ikB 2 lk w f (D\ lk t). 



(44) 



k=i 



In this sum, the first term with small eigenvalue Aio ~ 1 (in 2D) is explicitly separated 
from the remaining terms with large eigenvalues Aifc ~ ir 2 k 2 /£ 2 . As we mentioned in 
Sect. 3.2.2, the function Wf(p) slowly approaches 1 as p goes to 0, and has a power law 
decay (39) for p going to infinity. A large gap between Aio an d An creates a region of 
values p such that A^ 1 <C p <C A^ 1 , for which the first term in Eq. (44) is nearly constant, 
while the other terms can be neglected. This suggests an emergence of an intermediate 
diffusion regime with a nearly constant D(t), which is analogous to the tortuosity regime 
in porous media [97, 137]. This is a new, two-scale feature which could not be observed 
for one-scale domains such as the unit disk or sphere (Fig. 22). A complete analytical 
description developed in [128] allows one to study first the transition from the short-time 
diffusion to the intermediate (or tortuosity) regime at p around £ 2 , and then the transition 
to the ultimate long-time regime when p exceeds 1. 

In general, the emergence of an intermediate region with a constant D(t) would be 
a sign of a multi-scale geometry. Inversely, a significant separation in length scales in 
statistically homogeneous and isotropic porous media is expected to result in distinct 
constant regions in the behavior of D(t). Since these regions can in principle be observed 
by varying the diffusion time t, this observation suggests an experimental way for detecting 
multiple scales. However, a reliable interpretation of such measurements still requires a 
substantial study of the Laplace operator eigenbasis in porous structures. 

3.3 Residence times of reflected Brownian motion 

In the previous subsection, we focused on NMR-oriented applications, when the choice 
of the "observable" B(r) of Brownian motion to be a linear function of r was dictated 
by a widespread use of a linear magnetic field gradient in modern NMR scanners. The 
range of problems that can be tackled by matrix formalism is much broader. We consider 
here another important choice for the "observable", when B(r) is equal to the indicator 
function Ia (r) of a subset A of the confining domain: IU( r ) — 1 for r G A, and otherwise 
(Fig. 23). In this case, the random variable tp from Eq. (5) can be thought of as a "counter" 
which is turned on whenever the diffusing particles resides in A [138]. This is so-called 
residence or occupation time which is relevant for various diffusion-influenced reactions, 
for which the net outcome and the whole functioning of the system strongly depend on 
how long the diffusing particles remain in reactive zones (e.g., in a chemical reactor or 
biological cell) [139]. For instance, one can estimate the "trapping" time that particles 
spend in deep "fjord-like" pores of a catalyst. One can also mention surface relaxation 
processes in NMR [140] or kinetics of binding of ligands to a cell partially covered by 
receptors in microbiology [141]. 

The characteristic function E,{e iqtp } of the residence time (p is determined by Eq. (27), 
in which the matrix B is defined according to Eq. (20) with B(r) = !U(r): 



B, 



m,m 



i — 



/ 



dr u* m (r) I A (r) u m >(r) = 



dr u* m {r) u m ,{r). 



(45) 



n 



A 
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Figure 23: Diffusive motion of the particles restricted inside a bounded confining domain 
Q with reflecting boundary dQ. The residence time (p indicates how long a diffusing 
particle, started somewhere in the domain (e.g., in the subset A ), spends in a subset A 
until time t. The local time shows how long the particle spends in a close neighborhood 
(^-vicinity or e-sausage dQ £ ) of the boundary dQ. 



The sampling function p(r), determining the vector U, provides a way to weight or to 
condition the arrival points. In the simplest case with no conditioning, one uses p(r) = 1. 
If one aims to investigate the survival probability of only those particles that arrive in 
some subregion f2 of f2, one takes p(r) oc In ( r ) (Fig- 23). 

With the matrix formalism, one can easily assess much finer and more detailed statis- 
tics of residence times. For instance, to compute the residence time of the diffusing particle 
between two times < t\ < t 2 < t, a "counter" should be forced to turn on only during 
this period. Since the evolution of the system during two other periods [0, ti] and [t 2 ,t] 
is unperturbed (no "counting" or "interaction"), Eq. (27) can be modified as 

E{e iqtp } = (U e~ DAtl e-iDA-iqBWi-t!) e -DA(t-t 2 ) fjy 

Three matrix exponentials represent three successive time periods of the evolution. 
This "matrix product rule" can be applied in general, when one studies the residence 
time for a sequence of time intervals. Moreover, one can change the region of interest 
(set A) between different time periods in order to describe exchange processes between 
subsets (e.g., two pores of a medium). 

We employed the matrix representation of Sect. 3.1.2 to investigate the long-time be- 
havior of the moments of the residence time and to show that E{(p ra } is a polynomial 
of degree n as t goes to infinity [138]. The coefficients of this polynomial are expressed 
through the two governing matrices B and A, and vectors U and U, using a diagram- 
ing representation. In addition, we considered cumulant moments -< ip n >- defined as 
coefficients of the series 

lnE{e^} = J2^T * ^ y ■ 

n=l ^' 

We argued that all the cumulant moments behave as linear functions of t at long times: 

-< ip n ^~ b nA t + b nfi , 

where the coefficients b n ^ were explicitly given by using a diagrammic representation. 
This statement was explicitly demonstrated for the cumulant moments up to the order 
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Figure 24: The probability distribution of the local time of reflected Brownian motion on 
the unit interval at times t — 1,5, 10 (stars, triangles, and circles, respectively) obtained 
by Monte Carlo simulations, and its Gaussian approximation with mean b^i = 2 and 
variance 62,1 =2/3 (dashed-dotted, dashed, and solid lines, respectively). 



4 and checked numerically for several higher orders. Since the variance -< •p 2 y linearly 
increase in time, it is natural to renormalize p> by \/t. The first cumulant moment of 
the new random variable (p/\/i is still increasing in time, while its variance approaches a 
constant. In contrast, higher-order cumulant moments with n > 2 go to 0. In the long- 
time limit, the probability distribution of the normalized random variable <p/\/t becomes 
closer and closer to a Gaussian distribution with mean b\ t \\/t and variance 62,1, where 

= #0,0, &2,1 = 2 &Q,m^m- 

in 

As a consequence, the normalized residence time ip/y/i is getting closer and closer to a 
Gaussian variable as time t grows. These theoretical results are successfully confronted 
with Monte Carlo simulations (Fig. 24). 



3.3.1 Local time on the boundary 

The definition of reflected Brownian motion in Sect. 1.3 involves the local time £ t which 
characterizes the residence time on the boundary 

t 

£ t = \im- [ dsI d n £ (X s ), 

e^O e J 


where dfl £ is the e-vicinity of the boundary dQ: dQ £ = {r 6 Q : |r — dQ\ < e} 
(Fig. 23). The local time It is also related to the statistics of finite-distance reflections from 
the boundary, which is crucial for intermittent Brownian dynamics [70, 142-144]. The 
characteristic function of the local time is given by Eq. (27) with the matrix B determined 
by substituting B(r) = lonjz i n Eq. (45). When normalized by e, the integral over the 
e- vicinity dQ £ converges to the surface integral over the (smooth) boundary dQ, yielding 

&m,m> = J dr U *m( V ) «m'(r). (46) 

m 
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Here, the passage to the limit e — > 0, which is in general delicate and time-consuming 
for other numerical techniques is implemented intrinsically. This is one of the crucial 
advantages of the matrix formalism for computing local times. 

3.3.2 Surface relaxation mechanism 

The matrix formula (27) for the macroscopic signal exhibits the explicit dependence on 
the physical parameters D and 7/3, characterizing respectively diffusion and encoding. In 
turn, the dependence on the surface relaxivity W remains implicit, as being incorporated 
via the Robin boundary condition (9b) for the eigenfunctions. As a consequence, the 
eigenfunctions, the eigenvalues, and the governing matrices A and B have to be recalcu- 
lated for each new value of the surface relaxivity W. Since their computation is in general 
the most time-consuming step, such an implicit dependence on W can be considered as a 
drawback of matrix formalisms. 

This drawback can be overcome by introducing relaxation mechanisms via an inhomo- 
geneous distribution of relaxation rates B(r) in Eq. (4) with Neumann boundary condition. 
Surface relaxation can be implemented via a distribution B(r) localized near the bound- 
ary, e.g., B e (r) = £ _1 Ian e (r). In the limit e going to 0, the volume integral in Eq. (20), 
determining the governing matrix B, is reduced to the boundary integral (46). 

Since the mechanisms of gradient encoding and of bulk or surface relaxations are 
independent, their effects are simply superimposed as a linear combination of the cor- 
responding terms in Eq. (4). Consequently, the above expressions for the signal can be 
easily modified to include different attenuation mechanisms. For instance, Eq. (27) for 
the FID in the presence of surface relaxation becomes 

E = (U e -( DA + i 7/ 3e + w/ ^)t f j^ 

The advantage of this relation is the explicit dependence on all three physical parameters 
D, 7/5, and W. The structure of each term has a clear physical interpretation: DA de- 
scribes restricted diffusion, ry/3£> represents the dephasing, and WB S accounts for surface 
relaxation. One can similarly extend other matrix formulas for the CPMG sequence or 
any temporal profile. Nonuniform bulk or surface relaxations can also be incorporated. 

It is crucial to stress that here the eigenvalues and eigenfunctions are defined for 
Neumann boundary condition, whatever the value of surface relaxivity W is. As a con- 
sequence, these eigenfunctions, as well as the governing matrices A, B, and B s , depend 
only on the diffusion-confining geometry and have to be constructed only once for a given 
confining domain. The introduction of the matrix B s is therefore an alternative way 
to incorporate uniform surface relaxation. Moreover, it is a general frame for dealing 
with a nonuniform distribution B(r) of the relaxation rates, either in the bulk, or on the 
boundary. This concept is easily extendable for a superposition of various attenuation 
mechanisms. For instance, one can study the combined effect of the surface and bulk 
relaxations, gradient encoding, presence of dipolar magnetic field, etc. 
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Conclusion 



In this work, we presented theoretical and numerical results that we achieved in the 
course of the last five years. In order to study restricted diffusion in complex geometries, 
we combined probabilistic tools with spectral analysis. 

An implementation of Monte Carlo techniques that we specifically adapted to several 
model geometries allowed us to shed a new light onto diffusion in complex media. With 
these tools, we tackled such problems as multifractal properties of the harmonic measure in 
2D and 3D (Sect. 2.2), scaling properties of the spread harmonic measures and the role of 
the exploration length (Sect. 2.3), first passage statistics (Sect. 2.4), passivation processes 
(Sect. 2.5), diffusion-weighted imaging of the lungs (Sect. 2.6). These problems come 
from different fields (harmonic analysis, heterogeneous catalysis, heat transfer, nuclear 
magnetic resonance, physiology, medicine, etc.), in which geometrical complexity plays a 
central role. 

Monte Carlo techniques successfully answer the question how particles diffuse in a given 
medium. However, they may fail to explain why some features of diffusion are common 
while others are too specific or geometry-dependent. In other words, what differentiates 
various confining media in respect of diffusion? This question is tightly related to inverse 
and optimization problems which consist in determining or designing geometry according 
to some diffusion characteristics. A spectral approach provides a unifying mathematical 
"language" for tackling these problems. In fact, many diffusion characteristics can be 
expressed in terms of the Laplace operator eigenfunctions whose properties are there- 
fore a source of potentially important information. For instance, the interplay between 
these eigenfunctions and a linear magnetic field gradient yielded an intermediate diffusion 
regime for thin circular and spherical layers (Sect. 3.2.6). More generally, the spectral 
approach was successfully applied to retrieve, extend, and critically discuss various re- 
sults on restricted diffusion in NMR (Sect. 3.2). Residence time and other functionals 
of reflected Brownian motion could as well be analyzed (Sect. 3.3). The investigation 
of restricted diffusion in complex geometries based on the study of the properties of the 
underlying Laplace operator eigenbasis is our principal guideline. 

Many questions about diffusion in complex geometries remain open. On the one hand, 
a further study of restricted diffusion in model geometries (such as packs of spheres or 
cylinders) will reveal and help to better understand various features of this process. On 
the other hand, it is possible to apply the developed arsenal of theoretical and numeri- 
cal methods to complex structures in nature and industry: a three-dimensional porous 
morphology of a cement paste reconstructed by X-ray microtomography; skin interstitial 
space; the human placenta; the lungs, to name a few. Although a complete compre- 
hension of diffusion in natural geometries is irrealistic because of their complexity, this 
very complexity may help in averaging out specific features of the confining media and of 
transport processes. A practical goal of this study is to reveal relevant characteristics of 
the medium that would essentially determine and control diffusive transport. 
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